[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[gnuastro-commits] master 1a91475 2/2: Wrappers for GSL's 1D interpolati
|
From: |
Mohammad Akhlaghi |
|
Subject: |
[gnuastro-commits] master 1a91475 2/2: Wrappers for GSL's 1D interpolation added to library |
|
Date: |
Tue, 29 May 2018 21:25:06 -0400 (EDT) |
branch: master
commit 1a91475545152a77590411b5bab13491145b67f1
Author: Mohammad Akhlaghi <address@hidden>
Commit: Mohammad Akhlaghi <address@hidden>
Wrappers for GSL's 1D interpolation added to library
Two basic wrappers were added to make it easy for using GSL's 1D
interpolation functions with datasets stored in Gnuastro's `gal_data_t'.
Some minor corrections in the indexs of the manual were also found and are
also commited here.
Also, while working on text files, I noticed that the error message when
trying to read a text table wasn't too clear, so it is made more clear with
this commit.
---
NEWS | 6 +-
doc/gnuastro.texi | 251 +++++++++++++++++++++++++++++++++------
lib/gnuastro/interpolate.h | 20 +++-
lib/interpolate.c | 289 +++++++++++++++++++++++++++++++++++++++++++++
lib/txt.c | 2 +-
5 files changed, 525 insertions(+), 43 deletions(-)
diff --git a/NEWS b/NEWS
index 1f2ca39..db4fa6b 100644
--- a/NEWS
+++ b/NEWS
@@ -82,6 +82,8 @@ GNU Astronomy Utilities NEWS -*-
outline -*-
gal_eps_suffix_is_eps: Returns 1 if given suffix is EPS.
gal_eps_to_pt: Converts dataset size to PostScript points.
gal_eps_write: Writes a dataset into an EPS file.
+ gal_interpolate_1d_blank: Fill blank elements through interpolation.
+ gal_interpolate_1d_make_gsl_spline: Allocate and initalize `gsl_spline'.
gal_jpeg_name_is_jpeg: Returns 1 if given filename is JPEG.
gal_jpeg_suffix_is_jpeg: Returns 1 if given suffix is JPEG.
gal_jpeg_read: Reads input JPEG image into `gal_data_t'.
@@ -95,8 +97,8 @@ GNU Astronomy Utilities NEWS -*-
outline -*-
gal_pointer_allocate_mmap: Allocate space in a file, not in RAM.
gal_qsort_index_single_d: Sort indexs of single array in decreasing order.
gal_qsort_index_single_i: Sort indexs of single array in decreasing order.
- gal_qsort_index_multi_d: Sort indexs in multiple arrays (on different
threads).
- gal_qsort_index_multi_i: Sort indexs in multiple arrays (on different
threads).
+ gal_qsort_index_multi_d: Sort indexs in multiple arrays (different
threads).
+ gal_qsort_index_multi_i: Sort indexs in multiple arrays (different
threads).
gal_tiff_name_is_tiff: check if name contains a TIFF suffix.
gal_tiff_suffix_is_tiff: check if suffix is a TIFF suffix.
gal_tiff_dir_string_read: convert a string to a TIFF directory number.
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index d1f9924..e8f9ba5 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -1423,13 +1423,13 @@ tablets which allow you to do this.
@cindex Inconsistent results
According to Wikipedia ``a software bug is an error, flaw, failure, or
fault in a computer program or system that causes it to produce an
-incorrect or unexpected result, or to behave in unintended ways''. So
-when you see that a program is crashing, not reading your input
-correctly, giving the wrong results, or not writing your output
-correctly, you have found a bug. In such cases, it is best if you
-report the bug to the developers. The programs will also report bugs
-in known impossible situations (which are caused by something
-unexpected) and will ask the users to report the bug.
+incorrect or unexpected result, or to behave in unintended ways''. So when
+you see that a program is crashing, not reading your input correctly,
+giving the wrong results, or not writing your output correctly, you have
+found a bug. In such cases, it is best if you report the bug to the
+developers. The programs will also inform you if known impossible
+situations occur (which are caused by something unexpected) and will ask
+the users to report the bug issue.
@cindex Bug reporting
Prior to actually filing a bug report, it is best to search previous
@@ -4987,6 +4987,8 @@ line.
@item Debian-based OSs (@command{apt-get})
@cindex Debian
@cindex Ubuntu
address@hidden @command{apt-get}
address@hidden Advanced Packaging Tool (APT, 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
@@ -5002,6 +5004,9 @@ $ sudo apt-get install ghostscript libtool-bin
libjpeg-dev \
@item macOS (@command{brew})
@cindex macOS
address@hidden Homebrew
address@hidden MacPorts
address@hidden @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
@@ -5025,6 +5030,7 @@ $ brew install wcslib
@item Arch Linux (@command{pacman})
@cindex Arch Linux
address@hidden @command{pacman}
Arch Linux is a smaller GNU/Linux distribution. As discussed in
@url{https://en.wikipedia.org/wiki/Arch_Linux, Wikipedia}, it follows ``the
KISS principle ("keep it simple, stupid") as the general guideline, and
@@ -23677,7 +23683,7 @@ $ asttable --info table.fits
@end example
@end deftypefun
address@hidden {gal_data_t *} gal_table_read (char @code{*filename}, char
@code{*hdu}, gal_list_str_t @code{*cols}, int @code{searchin}, int
@code{ignorecase}, size_t @code{minmapsize}, size_t @code{colmatch})
address@hidden {gal_data_t *} gal_table_read (char @code{*filename}, char
@code{*hdu}, gal_list_str_t @code{*cols}, int @code{searchin}, int
@code{ignorecase}, size_t @code{minmapsize}, size_t @code{*colmatch})
Read the specified columns in a text file (named @code{filename}) into a
linked list of data structures. If the file is FITS, then @code{hdu} will
also be used, otherwise, @code{hdu} is ignored.
@@ -26872,7 +26878,7 @@ inside @code{indexs}. If @code{withrivers} is non-zero,
then pixels that
are immediately touching more than one positive value will be given a
@code{GAL_LABEL_RIVER} label.
address@hidden GNU C Library
address@hidden GNU C library
Note that the @code{indexs->array} is not re-allocated to its new size at
the address@hidden that according to the GNU C Library, even a
@code{realloc} to a smaller size can also cause a re-write of the whole
@@ -26942,34 +26948,35 @@ is much faster.
@node Interpolation, Git wrappers, Convolution functions, Gnuastro library
@subsection Interpolation (@file{interpolate.h})
-During data analysis, it often happens that parts of the data cannot be
-given a value. For example your image was saturated due to a very bright
-start and you have to mask that star's footprint. One other common
-situation in Gnuastro is when we do processing on tiles (for example to
-estimate the Sky value and its Standard deviation, see @ref{Sky
-value}). Some tiles must not be used for the estimation of the Sky value,
-for example because they cover a large galaxy. So we need to fill them in
-with blank values. But ultimately, we need a Sky value for every
-pixel. This job (assigning a value to blank element(s) based on their
-nearby neighbors with a value) is known as interpolation.
-
-There are many ways to do interpolation, but (mainly due to lack of time),
-currently Gnuastro only contains the (median of) nearest-neighbor
-method. The power of this method of interpolation is its non-parametric
-nature. The produced values are also always within the range of the known
-values and strong outliers do not get created. We will hopefully implement
-other methods too (wrappers around the GNU Scientific Library's very
-complete set of functions), but currently the developers are too busy. So
-if you do have the chance to help, your contribution would be very welcome
-and appreciated.
address@hidden Sky line
address@hidden Interpolation
+During data analysis, it happens that parts of the data cannot be given a
+value, but one is necessary for the higher-level analysis. For example a
+very bright star saturated part of your image and you need to fill in the
+saturated pixels with some values. Another common usage case are masked
+sky-lines in 1D specra that similarly need to be assigned a value for
+higher-level analysis. In other situations, you might want a value in an
+arbitrary point: between the elements/pixels where you have data. The
+functions described in this section are for such operations.
+
address@hidden GNU Scientific Library
+The parametric interpolations discussed below are wrappers around the
+interpolation functions of the GNU Scientific Library (or GSL, see @ref{GNU
+Scientific Library}). To identify the different GSL interpolation types,
+Gnuastro's @file{gnuastro/interpolate.h} header file contains macros that
+are discussed below. The GSL wrappers provided here are not yet complete
+because we are too busy. If you need them, please consider helping us in
+adding them to Gnuastro's library. Your would be very welcome and
+appreciated.
@deftypefun {gal_data_t *} gal_interpolate_close_neighbors (gal_data_t
@code{*input}, struct gal_tile_two_layer_params @code{*tl}, size_t
@code{numneighbors}, size_t @code{numthreads}, int @code{onlyblank}, int
@code{aslinkedlist})
Interpolate the values in the image using the median value of their
address@hidden closest neighbors. If @code{onlyblank} is non-zero,
-then only blank elements will be interpolated and pixels that already have
-a value will be left untouched. This function is multi-threaded and will
-run on @code{numthreads} threads (see @code{gal_threads_number} in
address@hidden programming}).
address@hidden closest neighbors. This function is non-parametric and
+thus agnostic to the input's number of dimension. If @code{onlyblank} is
+non-zero, then only blank elements will be interpolated and pixels that
+already have a value will be left untouched. This function is
+multi-threaded and will run on @code{numthreads} threads (see
address@hidden in @ref{Multithreaded programming}).
@code{tl} is Gnuastro's two later tessellation structure used to define
tiles over an image and is fully described in @ref{Tile grid}. When
@@ -26979,11 +26986,177 @@ properties, for example to not interpolate over
channel borders.
If several datasets have the same set of blank values, you don't need to
call this function multiple times. When @code{aslinkedlist} is non-zero,
-then @code{input} will be seen as a @ref{List of gal_data_t} and for all
-the same neighbor position checking will be done for all the datasets in
-the list. Of course, the values for each dataset will be different, so a
-different value will be written in the each dataset, but the neighbor
-checking that is the most CPU intensive part will only be done once.
+then @code{input} will be seen as a @ref{List of gal_data_t}. In this case,
+the same neighbors will be used for all the datasets in the list. Of
+course, the values for each dataset will be different, so a different value
+will be written in the each dataset, but the neighbor checking that is the
+most CPU intensive part will only be done once.
+
+This is a non-parametric and robust function for interpolation. The
+interpolated values are also always within the range of the non-blank
+values and strong outliers do not get created. However, this type of
+interpolation must be used with care when there are gradients. This is
+because it is non-parametric and if there aren't enough neighbors,
+step-like features can be created.
address@hidden deftypefun
+
address@hidden Macro GAL_INTERPOLATE_1D_INVALID
+This is just a place holder to manage errors.
address@hidden deffn
address@hidden Macro GAL_INTERPOLATE_1D_LINEAR
+[From GSL:] Linear interpolation. This interpolation method does not
+require any additional memory.
address@hidden deffn
address@hidden Macro GAL_INTERPOLATE_1D_POLYNOMIAL
address@hidden Polynomial interpolation
address@hidden Interpolation: Polynomial
+[From GSL:] Polynomial interpolation. This method should only be used for
+interpolating small numbers of points because polynomial interpolation
+introduces large oscillations, even for well-behaved datasets. The number
+of terms in the interpolating polynomial is equal to the number of points.
address@hidden deffn
address@hidden Macro GAL_INTERPOLATE_1D_CSPLINE
address@hidden Interpolation: Spline
address@hidden Cubic spline interpolation
address@hidden Spline (cubic) interpolation
+[From GSL:] Cubic spline with natural boundary conditions. The resulting
+curve is piecewise cubic on each interval, with matching first and second
+derivatives at the supplied data-points. The second derivative is chosen
+to be zero at the first point and last point.
address@hidden deffn
address@hidden Macro GAL_INTERPOLATE_1D_CSPLINE_PERIODIC
+[From GSL:] Cubic spline with periodic boundary conditions. The resulting
+curve is piecewise cubic on each interval, with matching first and second
+derivatives at the supplied data-points. The derivatives at the first and
+last points are also matched. Note that the last point in the data must
+have the same y-value as the first point, otherwise the resulting periodic
+interpolation will have a discontinuity at the boundary.
address@hidden deffn
address@hidden Macro GAL_INTERPOLATE_1D_AKIMA
address@hidden Interpolation: Akima spline
address@hidden Akima spline interpolation
address@hidden Spline (Akima) interpolation
+[From GSL:] Non-rounded Akima spline with natural boundary conditions. This
+method uses the non-rounded corner algorithm of Wodicka.
address@hidden deffn
address@hidden Macro GAL_INTERPOLATE_1D_AKIMA_PERIODIC
+[From GSL:] Non-rounded Akima spline with periodic boundary
+conditions. This method uses the non-rounded corner algorithm of Wodicka.
address@hidden deffn
address@hidden Macro GAL_INTERPOLATE_1D_STEFFEN
address@hidden Steffen interpolation
address@hidden Interpolation: Steffen
address@hidden Interpolation: monotonic
+[From GSL:] Steffen’s
address@hidden@url{http://adsabs.harvard.edu/abs/1990A%26A...239..443S}}
+guarantees the monotonicity of the interpolating function between the given
+data points. Therefore, minima and maxima can only occur exactly at the
+data points, and there can never be spurious oscillations between data
+points. The interpolated function is piecewise cubic in each interval. The
+resulting curve and its first derivative are guaranteed to be continuous,
+but the second derivative may be discontinuous.
address@hidden deffn
+
address@hidden {gsl_spline *} gal_interpolate_1d_make_gsl_spline (gal_data_t
@code{*X}, gal_data_t @code{*Y}, int @code{type_1d})
address@hidden GNU Scientific Library
+Allocate and initialize a GNU Scientific Library (GSL) 1D @code{gsl_spline}
+structure using the non-blank elements of @code{Y}. @code{type_1d}
+identifies the interpolation scheme and must be one of the
address@hidden macros defined above.
+
+If @code{X==NULL}, the X-axis is assumed to be integers starting from zero
+(the index of each element in @code{Y}). Otherwise, the values in @code{X}
+will be used to initialize the interpolation structure. Note that when
+given, @code{X} must @emph{not} contain any blank elements and it must be
+sorted (in increasing order).
+
+Each interpolation scheme needs a minimum number of elements to
+successfully operate. If the number of non-blank values in @code{Y} is less
+than this number, this function will return a @code{NULL} pointer.
+
+To be as generic and modular as possible, GSL's tools are low-level.
+Therefore before doing the interpolation, many steps are necessary (like
+preparing your dataset, then allocating and initializing
address@hidden). The metadata available in Gnuastro's @ref{Generic data
+container} make it easy to hide all those preparations within this
+function.
+
+Once @code{gsl_spline} has been initialized by this function, the
+interpolation can be evaluated for any X value within the non-blank range
+of the input using @code{gsl_spline_eval} or @code{gsl_spline_eval_e}.
+
+For example in the small program below, we read the first two columns of
+the table in @file{table.txt} and feed them to this function to later estimate
+the values in the second column for three selected points. You can use
address@hidden to compile and run this function, see @ref{Library demo
+programs} for more.
+
address@hidden first-in-first-out
address@hidden
+#include <stdio.h>
+#include <stdlib.h>
+#include <gnuastro/table.h>
+#include <gnuastro/interpolate.h>
+
+int
+main(void)
address@hidden
+ size_t i;
+ gal_data_t *X, *Y;
+ gsl_spline *spline;
+ gsl_interp_accel *acc;
+ gal_list_str_t *cols=NULL;
+
+ /* Change the values based on your input table. */
+ double address@hidden, 2.5, address@hidden;
+
+ /* Read the first two columns from `tab.txt'.
+ IMPORTANT: the list is first-in-first-out, so the output
+ column order is the inverse of the input order. */
+ gal_list_str_add(&cols, "1", 0);
+ gal_list_str_add(&cols, "2", 0);
+ Y=gal_table_read("table.txt", NULL, cols, GAL_TABLE_SEARCH_NAME,
+ 0, -1, NULL);
+ X=Y->next;
+
+ /* Allocate the GSL interpolation accelerator and make the
+ `gsl_spline' structure. */
+ acc=gsl_interp_accel_alloc();
+ spline=gal_interpolate_1d_make_gsl_spline(X, Y,
+ GAL_INTERPOLATE_1D_STEFFEN);
+
+ /* Calculate the respective value for all the given points,
+ if `spline' could be allocated. */
+ if(spline)
+ for(i=0; i<(sizeof points)/(sizeof *points); ++i)
+ printf("%f: %f\n", points[i],
+ gsl_spline_eval(spline, points[i], acc));
+
+ /* Clean up and return. */
+ gal_data_free(X);
+ gal_data_free(Y);
+ gsl_spline_free(spline);
+ gsl_interp_accel_free(acc);
+ gal_list_str_free(cols, 0);
+ return EXIT_SUCCESS;
address@hidden
address@hidden example
address@hidden deftypefun
+
address@hidden void gal_interpolate_1d_blank (gal_data_t @code{*in}, int
@code{type_1d})
+Fill the blank elements of @code{in} using the rest of the elements and the
+given interpolation. The interpolation scheme can be set through
address@hidden, which accepts any of the @code{GAL_INTERPOLATE_1D_*} macros
+above. The interpolation is internally done in 64-bit floating point type
+(@code{double}). However the evaluated/interpolated values (originally
+blank) will be written (in @code{in}) with its original numeric datatype,
+using C's standard type conversion.
+
+By definition, interpolation is only defined ``between'' valid
+points. Therefore, if any number of elements on the start or end of the 1D
+array are blank, those elements will not be interpolated and will remain
+blank. To see if any blank (non-interpolated) elements remain, you can use
address@hidden on @code{in} after this function is finished.
@end deftypefun
diff --git a/lib/gnuastro/interpolate.h b/lib/gnuastro/interpolate.h
index 0a5438d..53979a1 100644
--- a/lib/gnuastro/interpolate.h
+++ b/lib/gnuastro/interpolate.h
@@ -27,6 +27,7 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
must be included before the C++ preparations below */
#include <gnuastro/data.h>
#include <gnuastro/tile.h>
+#include <gsl/gsl_spline.h>
/* C++ Preparations */
#undef __BEGIN_C_DECLS
@@ -45,6 +46,20 @@ __BEGIN_C_DECLS /* From C++ preparations */
+/* Types of interpolation. */
+enum gal_interpolate_1D_types
+{
+ GAL_INTERPOLATE_1D_INVALID,
+
+ GAL_INTERPOLATE_1D_LINEAR,
+ GAL_INTERPOLATE_1D_POLYNOMIAL,
+ GAL_INTERPOLATE_1D_CSPLINE,
+ GAL_INTERPOLATE_1D_CSPLINE_PERIODIC,
+ GAL_INTERPOLATE_1D_AKIMA,
+ GAL_INTERPOLATE_1D_AKIMA_PERIODIC,
+ GAL_INTERPOLATE_1D_STEFFEN,
+};
+
gal_data_t *
@@ -53,8 +68,11 @@ gal_interpolate_close_neighbors(gal_data_t *input,
size_t numneighbors, size_t numthreads,
int onlyblank, int aslinkedlist);
+gsl_spline *
+gal_interpolate_1d_make_gsl_spline(gal_data_t *X, gal_data_t *Y, int type_1d);
-
+void
+gal_interpolate_1d_blank(gal_data_t *in, int type_1d);
__END_C_DECLS /* From C++ preparations */
diff --git a/lib/interpolate.c b/lib/interpolate.c
index 4d9dceb..29cfddc 100644
--- a/lib/interpolate.c
+++ b/lib/interpolate.c
@@ -45,6 +45,7 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
/*********************************************************************/
/******************** Nearest neighbor ********************/
+/*************** (Dimension agnostic) ****************/
/*********************************************************************/
/* These are bit-flags, so we're using hexadecimal notation. */
#define INTERPOLATE_FLAGS_NO 0
@@ -447,3 +448,291 @@ gal_interpolate_close_neighbors(gal_data_t *input,
gal_list_void_free(prm.ngb_vals, 1);
return prm.out;
}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+/*********************************************************************/
+/******************** 1D on grid ********************/
+/*********************************************************************/
+gsl_spline *
+gal_interpolate_1d_make_gsl_spline(gal_data_t *X, gal_data_t *Y, int type_1d)
+{
+ size_t i, c;
+ double *x, *y;
+ gal_data_t *Xd, *Yd;
+ gsl_spline *spline=NULL;
+ const gsl_interp_type *itype=NULL;
+ int Yhasblank=gal_blank_present(Y, 0);
+
+ /* A small sanity check. */
+ if(Y->ndim!=1)
+ error(EXIT_FAILURE, 0, "%s: input dataset is not 1D (it is %zuD)",
+ __func__, Y->ndim);
+ if(X)
+ {
+ if( gal_dimension_is_different(X, Y) )
+ error(EXIT_FAILURE, 0, "%s: when two inputs are given, they must "
+ "have the same dimensions. X has %zu elements, while Y has "
+ "%zu", __func__, X->size, Y->size);
+ if(gal_blank_present(X, 0))
+ error(EXIT_FAILURE, 0, "%s: the X dataset has blank elements",
+ __func__);
+ }
+
+ /* Set the interpolation type. */
+ switch(type_1d)
+ {
+ case GAL_INTERPOLATE_1D_LINEAR:
+ itype=gsl_interp_linear; break;
+ case GAL_INTERPOLATE_1D_POLYNOMIAL:
+ itype=gsl_interp_polynomial; break;
+ case GAL_INTERPOLATE_1D_CSPLINE:
+ itype=gsl_interp_cspline; break;
+ case GAL_INTERPOLATE_1D_CSPLINE_PERIODIC:
+ itype=gsl_interp_cspline_periodic; break;
+ case GAL_INTERPOLATE_1D_AKIMA:
+ itype=gsl_interp_akima; break;
+ case GAL_INTERPOLATE_1D_AKIMA_PERIODIC:
+ itype=gsl_interp_akima_periodic; break;
+ case GAL_INTERPOLATE_1D_STEFFEN:
+ itype=gsl_interp_steffen; break;
+ default:
+ error(EXIT_FAILURE, 0, "%s: code %d not recognizable for the GSL "
+ "interpolation type", __func__, type_1d);
+ }
+
+ /* Initializations. Note that if Y doesn't have any blank elements and is
+ already in `double' type, then we don't need to make a copy. */
+ Yd = ( (Yhasblank || Y->type!=GAL_TYPE_FLOAT64)
+ ? gal_data_copy_to_new_type(Y, GAL_TYPE_FLOAT64)
+ : Y );
+ Xd = ( X
+ /* Has to be `Yhasblank', we KNOW X doesn't have blank values. */
+ ? ( (Yhasblank || X->type!=GAL_TYPE_FLOAT64)
+ ? gal_data_copy_to_new_type(X, GAL_TYPE_FLOAT64)
+ : X )
+ : gal_data_alloc(NULL, GAL_TYPE_FLOAT64, 1, Y->dsize, NULL,
+ 0, -1, NULL, NULL, NULL) );
+
+ /* Fill in the X axis values while also removing NaN/blank elements. */
+ c=0;
+ x=Xd->array;
+ y=Yd->array;
+ for(i=0;i<Yd->size;++i)
+ if( !isnan(y[i]) )
+ {
+ y[ c ] = y[i];
+ x[ c++ ] = X ? x[i] : i;
+ }
+
+ /* Make sure we have enough valid points for interpolation. */
+ if( c>=gsl_interp_type_min_size(itype) )
+ {
+ spline=gsl_spline_alloc(itype, c);
+ gsl_spline_init(spline, x, y, c);
+ }
+ else
+ spline=NULL;
+
+ /* Clean up and return. */
+ if(Xd!=X) gal_data_free(Xd);
+ if(Yd!=Y) gal_data_free(Yd);
+ return spline;
+}
+
+
+
+
+
+/* Return 0 if all blanks were filled. */
+static int
+interpolate_1d_blank_write(gal_data_t *in, gsl_spline *spline,
+ gsl_interp_accel *acc)
+{
+ double tmp;
+ int hasblank;
+ uint8_t *su8 =in->array, *u8 =in->array, *u8f =u8 +in->size;
+ int8_t *si8 =in->array, *i8 =in->array, *i8f =i8 +in->size;
+ uint16_t *su16=in->array, *u16=in->array, *u16f=u16+in->size;
+ int16_t *si16=in->array, *i16=in->array, *i16f=i16+in->size;
+ uint32_t *su32=in->array, *u32=in->array, *u32f=u32+in->size;
+ int32_t *si32=in->array, *i32=in->array, *i32f=i32+in->size;
+ uint64_t *su64=in->array, *u64=in->array, *u64f=u64+in->size;
+ int64_t *si64=in->array, *i64=in->array, *i64f=i64+in->size;
+ float *sf32=in->array, *f32=in->array, *f32f=f32+in->size;
+ double *sf64=in->array, *f64=in->array, *f64f=f64+in->size;
+
+ switch(in->type)
+ {
+ case GAL_TYPE_UINT8:
+ do
+ if(*u8==GAL_BLANK_UINT8)
+ {
+ /* If the evaluation is good, this function will return 0. */
+ if( gsl_spline_eval_e(spline, u8-su8, acc, &tmp)==0 )
+ *u8=tmp;
+ else hasblank=1;
+ }
+ while(++u8<u8f);
+ break;
+ case GAL_TYPE_INT8:
+ do
+ if(*i8==GAL_BLANK_INT8)
+ {
+ if( gsl_spline_eval_e(spline, i8-si8, acc, &tmp)==0 )
+ *u16=tmp;
+ else hasblank=1;
+ }
+ while(++i8<i8f);
+ break;
+ case GAL_TYPE_UINT16:
+ do
+ if(*u16==GAL_BLANK_UINT16)
+ {
+ if( gsl_spline_eval_e(spline, u16-su16, acc, &tmp)==0 )
+ *u16=tmp;
+ else hasblank=1;
+ }
+ while(++u16<u16f);
+ break;
+ case GAL_TYPE_INT16:
+ do
+ if(*i16==GAL_BLANK_INT16)
+ {
+ if( gsl_spline_eval_e(spline, i16-si16, acc, &tmp)==0 )
+ *i16=tmp;
+ else hasblank=1;
+ }
+ while(++i16<i16f);
+ break;
+ case GAL_TYPE_UINT32:
+ do
+ if(*u32==GAL_BLANK_UINT32)
+ {
+ if( gsl_spline_eval_e(spline, u32-su32, acc, &tmp)==0 )
+ *u32=tmp;
+ else hasblank=1;
+ }
+ while(++u32<u32f);
+ break;
+ case GAL_TYPE_INT32:
+ do
+ if(*i32==GAL_BLANK_INT32)
+ {
+ if( gsl_spline_eval_e(spline, i32-si32, acc, &tmp)==0 )
+ *i32=tmp;
+ else hasblank=1;
+ }
+ while(++i32<i32f);
+ break;
+ case GAL_TYPE_UINT64:
+ do
+ if(*u64==GAL_BLANK_UINT64)
+ {
+ if( gsl_spline_eval_e(spline, u64-su64, acc, &tmp)==0 )
+ *u64=tmp;
+ else hasblank=1;
+ }
+ while(++u64<u64f);
+ break;
+ case GAL_TYPE_INT64:
+ do
+ if(*i64==GAL_BLANK_INT64)
+ {
+ if( gsl_spline_eval_e(spline, i64-si64, acc, &tmp)==0 )
+ *i64=tmp;
+ else hasblank=1;
+ }
+ while(++i64<i64f);
+ break;
+ case GAL_TYPE_FLOAT32:
+ do
+ if(isnan(*f32))
+ {
+ if( gsl_spline_eval_e(spline, f32-sf32, acc, &tmp)==0 )
+ *f32=tmp;
+ else hasblank=1;
+ }
+ while(++f32<f32f);
+ break;
+ case GAL_TYPE_FLOAT64:
+ do
+ if(isnan(*f64))
+ {
+ if( gsl_spline_eval_e(spline, f64-sf64, acc, f64) )
+ hasblank=1;
+ }
+ while(++f64<f64f);
+ break;
+ default:
+ error(EXIT_FAILURE, 0, "%s: code %d is not a recognized data type",
+ __func__, in->type);
+ }
+ return hasblank;
+}
+
+
+
+
+
+void
+gal_interpolate_1d_blank(gal_data_t *in, int type_1d)
+{
+ int hasblank;
+ gsl_spline *spline;
+ gsl_interp_accel *acc;
+
+ /* If there are no blank elements, just return. */
+ if(!gal_blank_present(in, 1)) return;
+
+ /* Initialize the necessary structures. */
+ spline=gal_interpolate_1d_make_gsl_spline(NULL, in, type_1d);
+
+ /* If any interpolation structure was actually made. */
+ if(spline)
+ {
+ /* Write the values in the blank elements. */
+ acc=gsl_interp_accel_alloc();
+ hasblank=interpolate_1d_blank_write(in, spline, acc);
+
+ /* For a check.
+ {
+ size_t i;
+ double *d;
+ gal_data_t *check=gal_data_copy_to_new_type(in, GAL_TYPE_FLOAT64);
+ d=check->array;
+ for(i=0;i<check->size;++i)
+ printf("%-10zu%f\n", i, d[i]);
+ gal_data_free(check);
+ }
+ */
+
+ /* Set the blank flags, note that `GAL_DATA_FLAG_BLANK_CH' is already set
+ by the top call to `gal_blank_present'. */
+ if(hasblank)
+ in->flag |= GAL_DATA_FLAG_HASBLANK;
+ else
+ in->flag &= ~GAL_DATA_FLAG_HASBLANK;
+
+ /* Clean up. */
+ gsl_spline_free(spline);
+ gsl_interp_accel_free(acc);
+ }
+}
diff --git a/lib/txt.c b/lib/txt.c
index 3005d53..097be54 100644
--- a/lib/txt.c
+++ b/lib/txt.c
@@ -546,7 +546,7 @@ txt_get_info(char *filename, int format, size_t *numdata,
size_t *dsize)
fp=fopen(filename, "r");
if(fp==NULL)
error(EXIT_FAILURE, errno, "%s: couldn't open to read as a plain "
- "text %s in %s", filename, format_err, __func__);
+ "text %s (from Gnuastro's `%s')", filename, format_err, __func__);
/* Allocate the space necessary to keep each line as we parse it. Note