[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[gnuastro-commits] master 213d1c3 106/125: gal_data_copy_* functions can
|
From: |
Mohammad Akhlaghi |
|
Subject: |
[gnuastro-commits] master 213d1c3 106/125: gal_data_copy_* functions can now work on tiles |
|
Date: |
Sun, 23 Apr 2017 22:36:48 -0400 (EDT) |
branch: master
commit 213d1c326ec52ceed5508ebfa2e23a9cc237f5f4
Author: Mohammad Akhlaghi <address@hidden>
Commit: Mohammad Akhlaghi <address@hidden>
gal_data_copy_* functions can now work on tiles
In order to continue work on the tiles and their interpolation some small
separate corrections were made. The most important is that the
`gal_tile_to_contiguous' function that was introduced in the previous
commit has now been removed. The `gal_data_copy_to_new_type' function (and
thus all the other `gal_data_copy...' functions can now work on tiles (copy
a tile to a contiguous patch of memory of any given type), so there was no
more need for `gal_tile_to_contiguous'.
The other changes are:
- The name `gal_wcs_array_from_wcsprm' was ambiguous so it has been
renamed to `gal_wcs_warp_matrix'. This was done as part of work to
prepare the WCS properties of a tile.
- `gal_data_copy_wcs' (in `gnuastro/data.h') has been replaced by
`gal_wcs_copy' in the `gnuastro/wcs.h' header.
- `gal_data_to_same_type' has been completely removed since it was very
specific to some special condition that doesn't exist any more (this
function isn't used any part of Gnuastro). The later defined
`gal_data_copy_to_new_type_free' (in combination with
`gal_data_copy_to_new_type') can be easily used instead.
- `gal_data_out_type' was just copy and pasted to a more relevant part of
`lib/data.c' (to do with type information). Until now it was near the
copy functions.
- `gal_wcs_on_tile' will now copy the WCS of the tile's allocated block
and correct it for use in the tile.
- The default settings for the starting and ending (inclusive) elements of
a tile were corrected (subtracted by one since it is meant to be
inclusive: to avoid bugs when changing dimensions which can lead to
large increments).
---
bin/statistics/statistics.c | 4 +-
bin/statistics/ui.c | 4 +
bin/warp/ui.c | 2 +-
lib/blank.c | 2 +-
lib/data.c | 275 ++++++++++++++++++--------------------------
lib/gnuastro/data.h | 13 +--
lib/gnuastro/tile.h | 8 +-
lib/gnuastro/wcs.h | 9 +-
lib/statistics.c | 2 +-
lib/tile.c | 53 ---------
lib/wcs.c | 93 ++++++++++++++-
11 files changed, 229 insertions(+), 236 deletions(-)
diff --git a/bin/statistics/statistics.c b/bin/statistics/statistics.c
index fa04ad4..c882f61 100644
--- a/bin/statistics/statistics.c
+++ b/bin/statistics/statistics.c
@@ -783,7 +783,7 @@ print_basics(struct statisticsparams *p)
to be found in place to save time/memory. But having a sorted array
can decrease the floating point accuracy of the standard deviation. So
we'll do the median calculation in the end.*/
- tmp=gal_statistics_mode(p->input, mirrdist, 0);
+ tmp=gal_statistics_mode(p->input, mirrdist, 1);
d=tmp->array;
if(d[2]>GAL_STATISTICS_MODE_GOOD_SYM)
{ /* Same format as `gal_data_write_to_string' */
@@ -793,7 +793,7 @@ print_basics(struct statisticsparams *p)
gal_data_free(tmp);
/* Find and print the median: */
- tmp=gal_statistics_median(p->input, 1);
+ tmp=gal_statistics_median(p->input, 0);
str=gal_data_write_to_string(tmp->array, tmp->type, 0);
printf(" %-*s %s\n", namewidth, "Median:", str);
gal_data_free(tmp);
diff --git a/bin/statistics/ui.c b/bin/statistics/ui.c
index 9498ab0..208d337 100644
--- a/bin/statistics/ui.c
+++ b/bin/statistics/ui.c
@@ -27,6 +27,7 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <error.h>
#include <stdio.h>
+#include <gnuastro/wcs.h>
#include <gnuastro/fits.h>
#include <gnuastro/tile.h>
#include <gnuastro/qsort.h>
@@ -806,6 +807,9 @@ ui_preparations(struct statisticsparams *p)
{
p->inputformat=INPUT_FORMAT_IMAGE;
p->input=gal_fits_img_read(p->inputname, cp->hdu, cp->minmapsize);
+ if(p->ontile)
+ gal_wcs_read(p->inputname, cp->hdu, 0, 0, &p->input->nwcs,
+ &p->input->wcs);
}
else
{
diff --git a/bin/warp/ui.c b/bin/warp/ui.c
index d9faf45..b958a0d 100644
--- a/bin/warp/ui.c
+++ b/bin/warp/ui.c
@@ -460,7 +460,7 @@ ui_matrix_make_align(struct warpparams *p, double *tmatrix)
/* Find the pixel scale along the two dimensions. Note that we will be
using the scale along the image X axis for both values. */
- w=gal_wcs_array_from_wcsprm(p->input->wcs);
+ w=gal_wcs_warp_matrix(p->input->wcs);
ps=gal_wcs_pixel_scale_deg(p->input->wcs);
diff --git a/lib/blank.c b/lib/blank.c
index 27e7d21..5281787 100644
--- a/lib/blank.c
+++ b/lib/blank.c
@@ -127,7 +127,7 @@ gal_blank_present(gal_data_t *input)
char **str, **strf;
size_t increment=0, num_increment=1;
gal_data_t *block=gal_tile_block(input);
- size_t start_end_inc[2]={0,block->size};
+ size_t start_end_inc[2]={0,block->size-1}; /* -1: this is INCLUSIVE. */
/* If there is nothing in the array (its size is zero), then return 0 (no
blank is present. */
diff --git a/lib/data.c b/lib/data.c
index 1ad6b21..a683ca0 100644
--- a/lib/data.c
+++ b/lib/data.c
@@ -33,7 +33,9 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <dirent.h>
#include <sys/mman.h>
+#include <gnuastro/wcs.h>
#include <gnuastro/data.h>
+#include <gnuastro/tile.h>
#include <gnuastro/blank.h>
#include <gnuastro/table.h>
#include <gnuastro/linkedlist.h>
@@ -248,6 +250,18 @@ gal_data_is_linked_list(uint8_t type)
+int
+gal_data_out_type(gal_data_t *first, gal_data_t *second)
+{
+ int ftype=gal_tile_block(first)->type;
+ int stype=gal_tile_block(second)->type;
+ return ftype > stype ? ftype : stype;
+}
+
+
+
+
+
@@ -400,36 +414,6 @@ gal_data_ptr_dist(void *earlier, void *later, uint8_t type)
-/* Copy the WCS structure from the input to the output structure. */
-void
-gal_data_copy_wcs(gal_data_t *in, gal_data_t *out)
-{
- if(in->wcs)
- {
- /* Allocate the output WCS structure. */
- errno=0;
- out->wcs=malloc(sizeof *out->wcs);
- if(out->wcs==NULL)
- error(EXIT_FAILURE, errno, "%zu bytes for out->wcs in "
- "gal_data_copy_wcs", sizeof *out->wcs);
-
- /* Initialize the allocated WCS structure. The WCSLIB manual says "On
- the first invokation, and only the first invokation, wcsprm::flag
- must be set to -1 to initialize memory management"*/
- out->wcs->flag=-1;
- wcsini(1, out->ndim, out->wcs);
-
- /* Copy the input WCS to the output WSC structure. */
- wcscopy(1, in->wcs, out->wcs);
- }
- else
- out->wcs=NULL;
-}
-
-
-
-
-
/* Allocate an array based on the value of type. Note that the argument
`size' is the number of elements, necessary in the array, the number of
bytes each element needs will be determined internaly by this function
@@ -634,12 +618,10 @@ gal_data_initialize(gal_data_t *data, void *array,
uint8_t type,
char *unit, char *comment)
{
size_t i;
- gal_data_t in;
/* Do the simple copying cases. For the display elements, set them all to
impossible (negative) values so if not explicitly set by later steps,
the default values are used if/when printing.*/
- data->wcs=NULL;
data->status=0;
data->next=NULL;
data->ndim=ndim;
@@ -653,10 +635,8 @@ gal_data_initialize(gal_data_t *data, void *array, uint8_t
type,
data->disp_fmt=data->disp_width=data->disp_precision=-1;
- /* Copy the WCS structure. Note that the `in' data structure was just
- defined to keep this pointer to call `gal_data_copy_wcs'. */
- in.wcs=wcs;
- gal_data_copy_wcs(&in, data);
+ /* Copy the WCS structure. */
+ data->wcs=gal_wcs_copy(wcs, ndim);
/* Allocate space for the dsize array, only if the data are to have any
@@ -1171,6 +1151,10 @@ data_copy_from_string(gal_data_t *from, gal_data_t *to)
if(from->type!=GAL_DATA_TYPE_STRING)
error(EXIT_FAILURE, 0, "`from' in `data_copy_from_string' must have "
"a string type.");
+ if(from->block)
+ error(EXIT_FAILURE, 0, "'data_copyt_from_string' doesn't currently "
+ "support tile inputs. Please contact us at %s so we can implement "
+ "this feature", PACKAGE_BUGREPORT);
/* Do the copying. */
for(i=0;i<from->size;++i)
@@ -1263,6 +1247,10 @@ data_copy_to_string(gal_data_t *from, gal_data_t *to)
if(to->type!=GAL_DATA_TYPE_STRING)
error(EXIT_FAILURE, 0, "`to' in `data_copy_to_string' must have a "
"string type");
+ if(from->block)
+ error(EXIT_FAILURE, 0, "'data_copyt_to_string' doesn't currently "
+ "support tile inputs. Please contact us at %s so we can implement "
+ "this feature", PACKAGE_BUGREPORT);
/* Do the copying */
switch(from->type)
@@ -1321,64 +1309,57 @@ data_copy_to_string(gal_data_t *from, gal_data_t *to)
-/* Copy to a new type for integers. */
-#define COPY_OTYPE_ITYPE_SET_INT(otype, itype) { \
- itype *ia=in->array, iblank; \
- otype *oa=out->array, *of=oa+out->size, oblank; \
- \
- /* Check if there are blank values in the input array and that */ \
- /* the types of the two structures are different. */ \
- if( in->type!=newtype && gal_blank_present(in) ) \
- { \
- /* Set the blank values */ \
- gal_blank_write(&iblank, in->type); \
- gal_blank_write(&oblank, newtype); \
+#define COPY_OT_IT_SET(OT, IT) { \
+ OT ob, *o=out->array; \
+ size_t increment=0, num_increment=1; \
+ IT ib, *ist, *i=in->array, *f=i+in->size; \
+ size_t mclen, contig_len=in->dsize[in->ndim-1]; \
+ size_t s_e_ind[2]={0,iblock->size-1}; /* -1: this is INCLUSIVE */ \
\
- /* Copy the input to the output. */ \
- do { *oa = *ia==iblank ? oblank : *ia; ia++; } while(++oa<of); \
- } \
+ /* If we are on a tile, the default values need to change. */ \
+ if(in!=iblock) \
+ ist=gal_tile_start_end_ind_inclusive(in, iblock, s_e_ind); \
\
- /* There were no blank elements in the input, or the input and */ \
- /* output have the same type. */ \
+ /* Constant preparations before the loop. */ \
+ if(iblock->type==out->type) \
+ mclen = in==iblock ? iblock->size : contig_len; \
else \
- do *oa=*ia++; while(++oa<of); \
- }
-
-
-
-
-
-/* Copy to a new type for floating point values. */
-#define COPY_OTYPE_ITYPE_SET_FLT(otype, itype) { \
- itype *ia=in->array, iblank; \
- otype *oa=out->array, *of=oa+out->size, oblank; \
+ { \
+ gal_blank_write(&ob, out->type); \
+ gal_blank_write(&ib, iblock->type); \
+ } \
\
- /* Check if there are blank values in the input array and that */ \
- /* the types of the two structures are different. */ \
- if( in->type!=newtype && gal_blank_present(in) ) \
+ /* Parse over the input and copy it. */ \
+ while( s_e_ind[0] + increment <= s_e_ind[1] ) \
{ \
- /* Set the blank values */ \
- gal_blank_write(&iblank, in->type); \
- gal_blank_write(&oblank, newtype); \
+ /* If we are on a tile, reset `i' and `f' for each round. */ \
+ if(in!=iblock) \
+ f = ( i = ist + increment ) + contig_len; \
\
- /* When the blank value isn't NaN, then we should use the */ \
- /* equal operator to check for blank values. */ \
- if( isnan(iblank) ) \
+ /* When the types are the same just use memcopy, otherwise, */ \
+ /* We'll have to read each number (and use internal */ \
+ /* conversion). */ \
+ if(iblock->type==out->type) \
{ \
- do { *oa = isnan(*ia) ? oblank : *ia; ia++; } \
- while(++oa<of); \
+ memcpy(o, i, mclen*gal_data_sizeof(iblock->type)); \
+ o += mclen; \
} \
else \
{ \
- do { *oa = *ia==iblank ? oblank : *ia; ia++; } \
- while(++oa<of); \
+ /* If the blank is a NaN value (only for floating point */ \
+ /* types), it will fail any comparison, so we'll exploit */ \
+ /* this property in such cases. For other cases, a */ \
+ /* `*i==ib' is enough. */ \
+ if(ib==ib) do *o++ = *i==ib ? ob : *i; while(++i<f); \
+ else do *o++ = *i!=*i ? ob : *i; while(++i<f); \
} \
- } \
\
- /* There were no blank elements in the input, or the input and */ \
- /* output have the same type. */ \
- else \
- do *oa=*ia++; while(++oa<of); \
+ /* Update the increment from the start of the input. */ \
+ increment += ( in==iblock ? iblock->size \
+ : gal_tile_block_increment(iblock, in->dsize, \
+ num_increment++, \
+ NULL) ); \
+ } \
}
@@ -1387,47 +1368,47 @@ data_copy_to_string(gal_data_t *from, gal_data_t *to)
/* gal_data_copy_to_new_type: Output type is set, now choose the input
type. */
-#define COPY_OTYPE_SET(otype) \
- switch(in->type) \
+#define COPY_OT_SET(OT) \
+ switch(iblock->type) \
{ \
case GAL_DATA_TYPE_UINT8: \
- COPY_OTYPE_ITYPE_SET_INT(otype, uint8_t); \
+ COPY_OT_IT_SET(OT, uint8_t); \
break; \
\
case GAL_DATA_TYPE_INT8: \
- COPY_OTYPE_ITYPE_SET_INT(otype, int8_t); \
+ COPY_OT_IT_SET(OT, int8_t); \
break; \
\
case GAL_DATA_TYPE_UINT16: \
- COPY_OTYPE_ITYPE_SET_INT(otype, uint16_t); \
+ COPY_OT_IT_SET(OT, uint16_t); \
break; \
\
case GAL_DATA_TYPE_INT16: \
- COPY_OTYPE_ITYPE_SET_INT(otype, int16_t); \
+ COPY_OT_IT_SET(OT, int16_t); \
break; \
\
case GAL_DATA_TYPE_UINT32: \
- COPY_OTYPE_ITYPE_SET_INT(otype, uint32_t); \
+ COPY_OT_IT_SET(OT, uint32_t); \
break; \
\
case GAL_DATA_TYPE_INT32: \
- COPY_OTYPE_ITYPE_SET_INT(otype, int32_t); \
+ COPY_OT_IT_SET(OT, int32_t); \
break; \
\
case GAL_DATA_TYPE_UINT64: \
- COPY_OTYPE_ITYPE_SET_INT(otype, uint64_t); \
+ COPY_OT_IT_SET(OT, uint64_t); \
break; \
\
case GAL_DATA_TYPE_INT64: \
- COPY_OTYPE_ITYPE_SET_INT(otype, int64_t); \
+ COPY_OT_IT_SET(OT, int64_t); \
break; \
\
case GAL_DATA_TYPE_FLOAT32: \
- COPY_OTYPE_ITYPE_SET_FLT(otype, float); \
+ COPY_OT_IT_SET(OT, float); \
break; \
\
case GAL_DATA_TYPE_FLOAT64: \
- COPY_OTYPE_ITYPE_SET_FLT(otype, double); \
+ COPY_OT_IT_SET(OT, double); \
break; \
\
case GAL_DATA_TYPE_STRING: \
@@ -1445,18 +1426,21 @@ data_copy_to_string(gal_data_t *from, gal_data_t *to)
\
default: \
error(EXIT_FAILURE, 0, "type code %d not recognized for " \
- "`in->type' in COPY_OTYPE_SET", in->type); \
+ "`in->type' in COPY_OT_SET", in->type); \
}
-/* Copy a given data structure to a new one (possibly with a new type). */
+/* Copy a given data structure to a new one with any type (for the
+ output). The input can be a tile, in which case the output will be a
+ contiguous patch of memory that has all the values within the input tile
+ in the requested type. */
gal_data_t *
gal_data_copy_to_new_type(gal_data_t *in, uint8_t newtype)
{
- gal_data_t *out;
+ gal_data_t *out, *iblock=gal_tile_block(in);
/* Allocate space for the output type */
out=gal_data_alloc(NULL, newtype, in->ndim, in->dsize, in->wcs,
@@ -1475,18 +1459,18 @@ gal_data_copy_to_new_type(gal_data_t *in, uint8_t
newtype)
*/
/* Fill in the output array: */
- switch(newtype)
+ switch(out->type)
{
- case GAL_DATA_TYPE_UINT8: COPY_OTYPE_SET(uint8_t); break;
- case GAL_DATA_TYPE_INT8: COPY_OTYPE_SET(int8_t); break;
- case GAL_DATA_TYPE_UINT16: COPY_OTYPE_SET(uint16_t); break;
- case GAL_DATA_TYPE_INT16: COPY_OTYPE_SET(int16_t); break;
- case GAL_DATA_TYPE_UINT32: COPY_OTYPE_SET(uint32_t); break;
- case GAL_DATA_TYPE_INT32: COPY_OTYPE_SET(int32_t); break;
- case GAL_DATA_TYPE_UINT64: COPY_OTYPE_SET(uint64_t); break;
- case GAL_DATA_TYPE_INT64: COPY_OTYPE_SET(int64_t); break;
- case GAL_DATA_TYPE_FLOAT32: COPY_OTYPE_SET(float); break;
- case GAL_DATA_TYPE_FLOAT64: COPY_OTYPE_SET(double); break;
+ case GAL_DATA_TYPE_UINT8: COPY_OT_SET( uint8_t ); break;
+ case GAL_DATA_TYPE_INT8: COPY_OT_SET( int8_t ); break;
+ case GAL_DATA_TYPE_UINT16: COPY_OT_SET( uint16_t ); break;
+ case GAL_DATA_TYPE_INT16: COPY_OT_SET( int16_t ); break;
+ case GAL_DATA_TYPE_UINT32: COPY_OT_SET( uint32_t ); break;
+ case GAL_DATA_TYPE_INT32: COPY_OT_SET( int32_t ); break;
+ case GAL_DATA_TYPE_UINT64: COPY_OT_SET( uint64_t ); break;
+ case GAL_DATA_TYPE_INT64: COPY_OT_SET( int64_t ); break;
+ case GAL_DATA_TYPE_FLOAT32: COPY_OT_SET( float ); break;
+ case GAL_DATA_TYPE_FLOAT64: COPY_OT_SET( double ); break;
case GAL_DATA_TYPE_STRING: data_copy_to_string(in, out); break;
case GAL_DATA_TYPE_BIT:
@@ -1495,12 +1479,12 @@ gal_data_copy_to_new_type(gal_data_t *in, uint8_t
newtype)
case GAL_DATA_TYPE_COMPLEX64:
error(EXIT_FAILURE, 0, "`gal_data_copy_to_new_type' currently doesn't "
"support copying to %s type",
- gal_data_type_as_string(newtype, 1));
+ gal_data_type_as_string(out->type, 1));
break;
default:
error(EXIT_FAILURE, 0, "type %d not recognized for "
- "for newtype in gal_data_copy_to_new_type", newtype);
+ "for `out->type' in gal_data_copy_to_new_type", out->type);
}
/* Return the created array */
@@ -1511,21 +1495,31 @@ gal_data_copy_to_new_type(gal_data_t *in, uint8_t
newtype)
+/* Copy the input data structure into a new type */
gal_data_t *
gal_data_copy_to_new_type_free(gal_data_t *in, uint8_t type)
{
- gal_data_t *out;
+ gal_data_t *out, *iblock=gal_tile_block(in);
/* In a general application, it might happen that the type is equal with
- the type of the input. Since the job of this function is to free the
- input data set, so and the user just wants one dataset after this
- function finishes, we can safely just return the input. */
- if(type==in->type)
+ the type of the input and the input isn't a tile. Since the job of
+ this function is to free the input dataset, and the user just wants
+ one dataset after this function finishes, we can safely just return
+ the input. */
+ if(type==iblock->type && in==iblock)
return in;
else
{
out=gal_data_copy_to_new_type(in, type);
- gal_data_free(in);
+ if(iblock==in)
+ gal_data_free(in);
+ else
+ fprintf(stderr, "#####\nWarning from "
+ "`gal_data_copy_to_new_type_free'\n#####\n The input "
+ "dataset is a tile, not a contiguous (fully allocated) "
+ "patch of memory. So it has not been freed. Please use "
+ "`gal_data_copy_to_new_type' to avoid this warning.\n"
+ "#####");
return out;
}
}
@@ -1534,53 +1528,14 @@ gal_data_copy_to_new_type_free(gal_data_t *in, uint8_t
type)
+/* Wrapper for `gal_data_copy_to_new_type', but will copy to the same type
+ as the input. Recall that if the input is a tile (a part of the input,
+ which is not-contiguous if it has more than one dimension), then the
+ output will have only the elements that cover the tile.*/
gal_data_t *
gal_data_copy(gal_data_t *in)
{
- return gal_data_copy_to_new_type(in, in->type);
-}
-
-
-
-
-
-int
-gal_data_out_type(gal_data_t *first, gal_data_t *second)
-{
- return first->type > second->type ? first->type : second->type;
-}
-
-
-
-
-
-/* The two input `f' and `s' datasets can be any type. But `of' and `os'
- will have type `type', if freeinputs is non-zero, then the input arrays
- will be freed if they needed to be changed to a new type. */
-void
-gal_data_to_same_type(gal_data_t *f, gal_data_t *s,
- gal_data_t **of, gal_data_t **os,
- uint8_t type, int freeinputs)
-{
- /* Change first dataset into the new type if necessary. */
- if( f->type != type )
- {
- *of=gal_data_copy_to_new_type(f, type);
- if(freeinputs)
- gal_data_free(f);
- }
- else
- *of=f;
-
- /* Change second dataset into the new type if necessary. */
- if( s->type != type )
- {
- *os=gal_data_copy_to_new_type(s, type);
- if(freeinputs)
- gal_data_free(s);
- }
- else
- *os=s;
+ return gal_data_copy_to_new_type(in, gal_tile_block(in)->type);
}
diff --git a/lib/gnuastro/data.h b/lib/gnuastro/data.h
index a7b8799..de9cc89 100644
--- a/lib/gnuastro/data.h
+++ b/lib/gnuastro/data.h
@@ -248,14 +248,14 @@ gal_data_type_max(uint8_t type, void *in);
int
gal_data_is_linked_list(uint8_t type);
+int
+gal_data_out_type(gal_data_t *first, gal_data_t *second);
+
/*********************************************************************/
/************* Size and allocation *******************/
/*********************************************************************/
-void
-gal_data_copy_wcs(gal_data_t *in, gal_data_t *out);
-
int
gal_data_dsize_is_different(gal_data_t *first, gal_data_t *second);
@@ -349,13 +349,6 @@ gal_data_copy_to_new_type_free(gal_data_t *in, uint8_t
type);
gal_data_t *
gal_data_copy(gal_data_t *in);
-int
-gal_data_out_type(gal_data_t *first, gal_data_t *second);
-
-void
-gal_data_to_same_type(gal_data_t *f, gal_data_t *s, gal_data_t **of,
- gal_data_t **os, uint8_t type, int freeinputs);
-
void
gal_data_copy_element_same_type(gal_data_t *input, size_t index, void *ptr);
diff --git a/lib/gnuastro/tile.h b/lib/gnuastro/tile.h
index 62a5f02..3f5e478 100644
--- a/lib/gnuastro/tile.h
+++ b/lib/gnuastro/tile.h
@@ -59,9 +59,6 @@ void *
gal_tile_start_end_ind_inclusive(gal_data_t *tile, gal_data_t *work,
size_t *start_end_inc);
-gal_data_t *
-gal_tile_to_contiguous(gal_data_t *input);
-
/***********************************************************************/
@@ -157,12 +154,13 @@ gal_tile_function_on_threads(gal_data_t *tiles, void
*(*meshfunc)(void *),
See `lib/statistics.c' for some example applications of this function.
*/
#define GAL_TILE_PARSE_OPERATE(IT, OT, OP, IN, OUT, PARSE_OUT) { \
+ size_t increment=0, num_increment=1; \
OT *ost, *o = OUT ? OUT->array : NULL; \
gal_data_t *iblock=gal_tile_block(IN); \
IT b, *st, *i=IN->array, *f=i+IN->size; \
gal_data_t *oblock = OUT ? gal_tile_block(OUT) : NULL; \
- size_t increment=0, num_increment=1, s_e_ind[2]={0,iblock->size}; \
int hasblank=gal_blank_present(IN), parse_out=(OUT && PARSE_OUT); \
+ size_t s_e_ind[2]={0,iblock->size-1}; /* -1: this is INCLUSIVE */ \
\
/* Write the blank value for the input type into `b'). Note that */ \
/* a tile doesn't necessarily have to have a type. */ \
@@ -187,7 +185,7 @@ gal_tile_function_on_threads(gal_data_t *tiles, void
*(*meshfunc)(void *),
/* Go over contiguous patches of memory. */ \
while( s_e_ind[0] + increment <= s_e_ind[1] ) \
{ \
- /* If we are on a tile, reset `a' and `af'. Also, if there */ \
+ /* If we are on a tile, reset `i' and `f'. Also, if there */ \
/* is more than one element in `OUT', then set that. Note */ \
/* that it is upto the caller to increment `o' in `OP'. */ \
if(IN!=iblock) \
diff --git a/lib/gnuastro/wcs.h b/lib/gnuastro/wcs.h
index 117df70..d82c93e 100644
--- a/lib/gnuastro/wcs.h
+++ b/lib/gnuastro/wcs.h
@@ -28,6 +28,7 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <fitsio.h>
#include <wcslib/wcs.h>
+#include <gnuastro/data.h>
/* C++ Preparations */
@@ -69,8 +70,14 @@ gal_wcs_read(char *filename, char *hdu, size_t hstartwcs,
/**************************************************************/
/********** Utilities ************/
/**************************************************************/
+struct wcsprm *
+gal_wcs_copy(struct wcsprm *in, size_t ndim);
+
+void
+gal_wcs_on_tile(gal_data_t *tile);
+
double *
-gal_wcs_array_from_wcsprm(struct wcsprm *wcs);
+gal_wcs_warp_matrix(struct wcsprm *wcs);
void
gal_wcs_decompose_pc_cdelt(struct wcsprm *wcs);
diff --git a/lib/statistics.c b/lib/statistics.c
index 5669cc9..db4a3af 100644
--- a/lib/statistics.c
+++ b/lib/statistics.c
@@ -1293,7 +1293,7 @@ gal_statistics_no_blank_sorted(gal_data_t *input, int
inplace)
if(input->block)
{
/* Copy the input into a contiguous patch of memory. */
- contig=gal_tile_to_contiguous(input);
+ contig=gal_data_copy(input);
/* When the data was a tile, we have already copied the array into a
separate allocated space. So to avoid any further copying, we will
diff --git a/lib/tile.c b/lib/tile.c
index a40ef0b..2981f0a 100644
--- a/lib/tile.c
+++ b/lib/tile.c
@@ -189,59 +189,6 @@ gal_tile_start_end_ind_inclusive(gal_data_t *tile,
gal_data_t *work,
-/* Return a contiguous patch of memory with the same contents as the
- tile. If the input is not actually a tile, this function will return the
- actual input untouched. */
-#define TO_CONTIGUOUS(IT) { \
- IT *i, *o=out->array, *f, *st; \
- st=gal_tile_start_end_ind_inclusive(input, block, s_e_ind); \
- while( s_e_ind[0] + increment <= s_e_ind[1] ) \
- { \
- f = ( i = st + increment ) + input->dsize[input->ndim-1]; \
- do *o++=*i++; while(i<f); \
- increment += gal_tile_block_increment(block, input->dsize, \
- num_increment++, NULL); \
- } \
- }
-gal_data_t *
-gal_tile_to_contiguous(gal_data_t *input)
-{
- gal_data_t *out, *block=gal_tile_block(input);
- size_t s_e_ind[2], increment=0, num_increment=1;
-
- /* Check if this is actually a tile. */
- if(input->block)
- {
- /* Allocate the contiguous block of memory. */
- out=gal_data_alloc(NULL, block->type, input->ndim, input->dsize,
- NULL, 0, input->minmapsize, NULL, input->unit,
- NULL);
-
- /* Copy the tile's contents to the contiguous patch of memory. */
- switch(block->type)
- {
- case GAL_DATA_TYPE_UINT8: TO_CONTIGUOUS( uint8_t ); break;
- case GAL_DATA_TYPE_INT8: TO_CONTIGUOUS( int8_t ); break;
- case GAL_DATA_TYPE_UINT16: TO_CONTIGUOUS( uint16_t ); break;
- case GAL_DATA_TYPE_INT16: TO_CONTIGUOUS( int16_t ); break;
- case GAL_DATA_TYPE_UINT32: TO_CONTIGUOUS( uint32_t ); break;
- case GAL_DATA_TYPE_INT32: TO_CONTIGUOUS( int32_t ); break;
- case GAL_DATA_TYPE_UINT64: TO_CONTIGUOUS( uint64_t ); break;
- case GAL_DATA_TYPE_INT64: TO_CONTIGUOUS( int64_t ); break;
- case GAL_DATA_TYPE_FLOAT32: TO_CONTIGUOUS( float ); break;
- case GAL_DATA_TYPE_FLOAT64: TO_CONTIGUOUS( double ); break;
- }
- }
- else out=input;
-
- /* Return. */
- return out;
-}
-
-
-
-
-
diff --git a/lib/wcs.c b/lib/wcs.c
index e881110..8957f2e 100644
--- a/lib/wcs.c
+++ b/lib/wcs.c
@@ -34,7 +34,9 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <gsl/gsl_linalg.h>
#include <gnuastro/wcs.h>
+#include <gnuastro/tile.h>
#include <gnuastro/fits.h>
+#include <gnuastro/multidim.h>
@@ -184,8 +186,95 @@ gal_wcs_read(char *filename, char *hdu, size_t hstartwcs,
/**************************************************************/
/********** Utilities ************/
/**************************************************************/
+/* Copy a given WSC structure into another one. */
+struct wcsprm *
+gal_wcs_copy(struct wcsprm *in, size_t ndim)
+{
+ struct wcsprm *out;
+
+ /* If the input WCS is NULL, return a NULL WCS. */
+ if(in)
+ {
+ /* Allocate the output WCS structure. */
+ errno=0;
+ out=malloc(sizeof *out);
+ if(out==NULL)
+ error(EXIT_FAILURE, errno, "%zu bytes for out in `gal_wcs_copy'",
+ sizeof *out);
+
+ /* Initialize the allocated WCS structure. The WCSLIB manual says "On
+ the first invokation, and only the first invokation, wcsprm::flag
+ must be set to -1 to initialize memory management"*/
+ out->flag=-1;
+ wcsini(1, ndim, out);
+
+ /* Copy the input WCS to the output WSC structure. */
+ wcscopy(1, in, out);
+ }
+ else
+ out=NULL;
+
+ /* Return the final output. */
+ return out;
+}
+
+
+
+
+
+/* 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
+ creating the tiles they don't have any WCS structure. When needed, this
+ function can be used to add a WCS structure to the tile by copying the
+ WCS structure of its block and correcting its starting points. If the
+ tile already has a WCS structure, this function won't do anything.*/
+void
+gal_wcs_on_tile(gal_data_t *tile)
+{
+ size_t i, start_ind, ndim=tile->ndim;
+ gal_data_t *block=gal_tile_block(tile);
+ size_t *coord=gal_data_malloc_array(GAL_DATA_TYPE_SIZE_T, ndim);
+
+ /* If the tile already has a WCS structure, don't do anything. */
+ if(tile->wcs) return;
+ else
+ {
+ /* Copy the block's WCS into the tile. */
+ tile->wcs=gal_wcs_copy(block->wcs, block->ndim);
+
+ /* 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);
+
+ /* Correct the copied WCS structure. Note that crpix is indexed in
+ the FITS/Fortran order while coord is ordered in C, it also starts
+ counting from 1, not zero. */
+ for(i=0;i<ndim;++i)
+ tile->wcs->crpix[i] -= coord[ndim-1-i];
+ /*
+ printf("start_ind: %zu\n", start_ind);
+ printf("coord: %zu, %zu\n", coord[1]+1, coord[0]+1);
+ printf("CRPIX: %f, %f\n", tile->wcs->crpix[0], tile->wcs->crpix[1]);
+ */
+ }
+
+ /* Clean up. */
+ free(coord);
+}
+
+
+
+
+
+/* Return the Warping matrix of the given WCS structure. This will be the
+ final matrix irrespective of the type of storage in the WCS
+ structure. Recall that the FITS standard has several methods to store
+ the matrix, which is up to this function to account for and return the
+ final matrix. The output is an allocated DxD matrix where `D' is the
+ number of dimensions. */
double *
-gal_wcs_array_from_wcsprm(struct wcsprm *wcs)
+gal_wcs_warp_matrix(struct wcsprm *wcs)
{
double *out;
size_t i, j, size=wcs->naxis*wcs->naxis;
@@ -333,7 +422,7 @@ gal_wcs_pixel_scale_deg(struct wcsprm *wcs)
/* Write the full matrix into an array, irrespective of what type it was
stored in the wcsprm structure (`PCi_j' style or `CDi_j' style). */
- a=gal_wcs_array_from_wcsprm(wcs);
+ a=gal_wcs_warp_matrix(wcs);
/* Fill in the necessary GSL vector and Matrix structures. */
S.size=n; S.stride=1; S.data=pixscale;
- [gnuastro-commits] master 31b3b76 099/125: Spatial convolution now using new tessellation, (continued)
- [gnuastro-commits] master 31b3b76 099/125: Spatial convolution now using new tessellation, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master c0dd6db 102/125: The noedgecorrection option in spatial domain convolution, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 0dfdc34 069/125: Option values directly set in library, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 95ba655 100/125: Spatial convolution with new tessellation, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 0dec9ac 110/125: Allow for floating point errors in Crop comparisons, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master db4c16d 076/125: Pulled in all the recent work in master, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 64f8db5 080/125: ConvertType working except for text inputs and outputs, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 4d853b4 112/125: Old neighbors.h header removed, used new in MakeProfiles, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 9959e45 004/125: Arithmetic using new data structure, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 0a2fd0e 089/125: New implementations for histogram and CFP, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 213d1c3 106/125: gal_data_copy_* functions can now work on tiles,
Mohammad Akhlaghi <=
- [gnuastro-commits] master 4f6a8c4 095/125: Fits prints HDU info of FITS file with no options, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 63043e5 105/125: Statistical functions can accept tiled input, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 4839a11 104/125: Preparations for tile interpolation, Mohammad Akhlaghi, 2017/04/23
- [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