casacore
InterpolateArray1D.h
Go to the documentation of this file.
1 //# Interpolate1DArray.h: Interpolation in last dimension of an Array
2 //# Copyright (C) 1997,1999,2000,2001
3 //# Associated Universities, Inc. Washington DC, USA.
4 //#
5 //# This library is free software; you can redistribute it and/or modify it
6 //# under the terms of the GNU Library General Public License as published by
7 //# the Free Software Foundation; either version 2 of the License, or (at your
8 //# option) any later version.
9 //#
10 //# This library is distributed in the hope that it will be useful, but WITHOUT
11 //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
13 //# License for more details.
14 //#
15 //# You should have received a copy of the GNU Library General Public License
16 //# along with this library; if not, write to the Free Software Foundation,
17 //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
18 //#
19 //# Correspondence concerning AIPS++ should be addressed as follows:
20 //# Internet email: aips2-request@nrao.edu.
21 //# Postal address: AIPS++ Project Office
22 //# National Radio Astronomy Observatory
23 //# 520 Edgemont Road
24 //# Charlottesville, VA 22903-2475 USA
25 //#
26 //# $Id: Interpolate1DArray.h,v 8.1 1997/05/21 22:59:29 rm
27 
28 #ifndef SCIMATH_INTERPOLATEARRAY1D_H
29 #define SCIMATH_INTERPOLATEARRAY1D_H
30 
31 #include <casacore/casa/aips.h>
32 #include <casacore/casa/Arrays/ArrayFwd.h>
33 
34 namespace casacore { //# NAMESPACE CASACORE - BEGIN
35 
36 template <class T> class PtrBlock;
37 template <class T> class Block;
38 
39 // <summary> Interpolate in one dimension </summary>
40 
41 // <use visibility=export>
42 
43 // <reviewed reviewer="" date="" tests="" demos="">
44 // </reviewed>
45 
46 // <prerequisite>
47 // <li> <linkto class=Array>Array</linkto>
48 // <li> <linkto class=Vector>Vector</linkto>
49 // </prerequisite>
50 
51 // <etymology>
52 // The InterpolateArray1D class does interpolation in one dimension of
53 // an Array only.
54 // </etymology>
55 
56 // <synopsis>
57 // This class will, given the abscissa and ordinates of a set of one
58 // dimensional data, interpolate on this data set giving the value at any
59 // specified ordinate. It will extrapolate if necessary, but this is will
60 // usually give a poor result. There is no requirement for the ordinates to
61 // be regularly spaced, however they do need to be sorted and each
62 // abscissa should have a unique value.
63 //
64 // Interpolation can be done using the following methods:
65 // <ul>
66 // <li> Nearest Neighbour
67 // <li> Linear (default unless there is only one data point)
68 // <li> Cubic Polynomial
69 // <li> Natural Cubic Spline
70 // </ul>
71 //
72 // The abscissa must be a simple type (scalar value) that
73 // can be ordered. ie. an uInt, Int, Float or Double (not Complex). The
74 // ordinate can be an Array of any data type that has addition, and
75 // subtraction defined as well as multiplication by a scalar of the abcissa
76 // type.
77 // So the ordinate can be complex numbers, where the interpolation is done
78 // separately on the real and imaginary components.
79 // Use of Arrays as the the Range type is discouraged, operations will
80 // be very slow, it would be better to construct a single higher dimensional
81 // array that contains all the data.
82 //
83 // Note: this class (and these docs) are heavily based on the
84 // <linkto class=Interpolate1D>Interpolate1D</linkto>
85 // class in aips/Functionals. That class proved to be
86 // too slow for interpolation of large data volumes (i.e. spectral line
87 // visibility datasets) mainly due to the interface which forced the
88 // creation of large numbers of temporary Vectors and Arrays.
89 // This class is 5-10 times faster than Interpolate1D in cases where
90 // large amounts of data are to be interpolated.
91 // </synopsis>
92 
93 // <example>
94 // This code fragment does cubic interpolation on (xin,yin) pairs to
95 // produce (xout,yout) pairs.
96 // <srcblock>
97 // Vector<Float> xin(4); indgen(xin);
98 // Vector<Double> yin(4); indgen(yin); yin = yin*yin*yin;
99 // Vector<Float> xout(20);
100 // for (Int i=0; i<20; i++) xout(i) = 1 + i*0.1;
101 // Vector<Double> yout;
102 // InterpolateArray1D<Float, Double>::interpolate(yout, xout, xin, yin,
103 // InterpolateArray1D<Float,Double>::cubic);
104 // </srcblock>
105 // </example>
106 
107 // <motivation>
108 // This class was motivated by the need to interpolate visibilities
109 // in frequency to allow selection and gridding in velocity space
110 // with on-the-fly doppler correction.
111 // </motivation>
112 
113 // <templating arg=Domain>
114 // <li> The Domain class must be a type that can be ordered in a mathematical
115 // sense. This includes uInt, Int, Float, Double, but not Complex.
116 // </templating>
117 
118 // <templating arg=Range>
119 // <li> The Range class must have addition and subtraction of Range objects with
120 // each other as well as multiplication by a scalar defined. Besides the
121 // scalar types listed above this includes Complex, DComplex, and Arrays of
122 // any of these types. Use of Arrays is discouraged however.
123 // </templating>
124 
125 // <thrown>
126 // <li> AipsError
127 // </thrown>
128 // <todo asof="1997/06/17">
129 // <li> Implement flagging in cubic and spline interpolation
130 // </todo>
131 
132 
133 template <class Domain, class Range>
135 {
136 public:
137  // Interpolation methods
139  // nearest neighbour
141  // linear
143  // cubic
145  // cubic spline
146  spline
147  };
148 
149  // Interpolate in the last dimension of array yin whose x coordinates
150  // along this dimension are given by xin.
151  // Output array yout has interpolated values for x coordinates xout.
152  // E.g., interpolate a Cube(pol,chan,time) in the time direction, all
153  // values in the pol-chan plane are interpolated to produce the output
154  // pol-chan plane.
155  static void interpolate(Array<Range>& yout,
156  const Vector<Domain>& xout,
157  const Vector<Domain>& xin,
158  const Array<Range>& yin,
159  Int method);
160 
161  // deprecated version of previous function using Blocks - no longer needed
162  // now that Vector has a fast index operator [].
163  static void interpolate(Array<Range>& yout,
164  const Block<Domain>& xout,
165  const Block<Domain>& xin,
166  const Array<Range>& yin,
167  Int method);
168 
169  // Interpolate in the last dimension of array yin whose x coordinates
170  // along this dimension are given by xin.
171  // Output array yout has interpolated values for x coordinates xout.
172  // This version handles flagged data in a simple way: all outputs
173  // depending on a flagged input are flagged.
174  // If goodIsTrue==True, then that means
175  // a good data point has a flag value of True (usually for
176  // visibilities, good is False and for images good is True)
177  // If extrapolate==False, then xout points outside the range of xin
178  // will always be marked as flagged.
179  // TODO: implement flags for cubic and spline (presently input flags
180  // are copied to output).
181  static void interpolate(Array<Range>& yout,
182  Array<Bool>& youtFlags,
183  const Vector<Domain>& xout,
184  const Vector<Domain>& xin,
185  const Array<Range>& yin,
186  const Array<Bool>& yinFlags,
187  Int method,
188  Bool goodIsTrue=False,
189  Bool extrapolate=False);
190 
191  // deprecated version of previous function using Blocks - no longer needed
192  // now that Vector has a fast index operator [].
193  static void interpolate(Array<Range>& yout,
194  Array<Bool>& youtFlags,
195  const Block<Domain>& xout,
196  const Block<Domain>& xin,
197  const Array<Range>& yin,
198  const Array<Bool>& yinFlags,
199  Int method,
200  Bool goodIsTrue=False,
201  Bool extrapolate=False);
202 
203  // Interpolate in the middle axis in 3D array (yin) whose x coordinates along the
204  // this dimension are given by xin.
205  // Interpolate a Cube(pol,chan,time) in the chan direction.
206  // Currently only linear interpolation method is implemented.
207  // TODO: add support for nearest neiborhood, cubic, and cubic spline.
208  static void interpolatey(Cube<Range>& yout,
209  const Vector<Domain>& xout,
210  const Vector<Domain>& xin,
211  const Cube<Range>& yin,
212  Int method);
213 
214  // Interpolate in the middle dimension of 3D array yin whose x coordinates
215  // along this dimension are given by xin.
216  // Output array yout has interpolated values for x coordinates xout.
217  // This version handles flagged data in a simple way: all outputs
218  // depending on a flagged input are flagged.
219  // If goodIsTrue==True, then that means
220  // a good data point has a flag value of True (usually for
221  // visibilities, good is False and for images good is True)
222  // If extrapolate==False, then xout points outside the range of xin
223  // will always be marked as flagged.
224  // Currently only linear interpolation method is implemented.
225  // TODO: add support for nearest neiborhood, cubic, and cubic spline.
226  static void interpolatey(Cube<Range>& yout,
227  Cube<Bool>& youtFlags,
228  const Vector<Domain>& xout,
229  const Vector<Domain>& xin,
230  const Cube<Range>& yin,
231  const Cube<Bool>& yinFlags,
232  Int method,
233  Bool goodIsTrue=False,
234  Bool extrapolate=False);
235 
236 private:
237  // Interpolate the y-vectors of length ny from x values xin to xout.
238  static void interpolatePtr(PtrBlock<Range*>& yout,
239  Int ny,
240  const Vector<Domain>& xout,
241  const Vector<Domain>& xin,
242  const PtrBlock<const Range*>& yin,
243  Int method);
244 
245  // Interpolate the y-vectors of length ny from x values xin to xout.
246  // Take flagging into account
247  static void interpolatePtr(PtrBlock<Range*>& yout,
248  PtrBlock<Bool*>& youtFlags,
249  Int ny,
250  const Vector<Domain>& xout,
251  const Vector<Domain>& xin,
252  const PtrBlock<const Range*>& yin,
253  const PtrBlock<const Bool*>& yinFlags,
254  Int method, Bool goodIsTrue,
255  Bool extrapolate);
256 
257  // Interpolate along yaxis
259  Int na,
260  Int nb,
261  Int nc,
262  const Vector<Domain>& xout,
263  const Vector<Domain>& xin,
264  const PtrBlock<const Range*>& yin,
265  Int method);
266 
267  // Take flagging into account
268  static void interpolateyPtr(PtrBlock<Range*>& yout,
269  PtrBlock<Bool*>& youtFlags,
270  Int na,
271  Int nb,
272  Int nc,
273  const Vector<Domain>& xout,
274  const Vector<Domain>& xin,
275  const PtrBlock<const Range*>& yin,
276  const PtrBlock<const Bool*>& yinFlags,
277  Int method, Bool goodIsTrue,
278  Bool extrapolate);
279 
280  // Interpolate the y-vectors of length ny from x values xin to xout
281  // using polynomial interpolation with specified order.
283  Int ny,
284  const Vector<Domain>& xout,
285  const Vector<Domain>& xin,
286  const PtrBlock<const Range*>& yin,
287  Int order);
288 
289 };
290 
291 
292 
293 } //# NAMESPACE CASACORE - END
294 
295 #ifndef CASACORE_NO_AUTO_TEMPLATES
296 #include <casacore/scimath/Mathematics/InterpolateArray1D.tcc>
297 #endif //# CASACORE_NO_AUTO_TEMPLATES
298 #endif
InterpolationMethod
Interpolation methods.
@ nearestNeighbour
nearest neighbour
static void interpolateyPtr(PtrBlock< Range * > &yout, PtrBlock< Bool * > &youtFlags, Int na, Int nb, Int nc, const Vector< Domain > &xout, const Vector< Domain > &xin, const PtrBlock< const Range * > &yin, const PtrBlock< const Bool * > &yinFlags, Int method, Bool goodIsTrue, Bool extrapolate)
Take flagging into account.
static void interpolatey(Cube< Range > &yout, Cube< Bool > &youtFlags, const Vector< Domain > &xout, const Vector< Domain > &xin, const Cube< Range > &yin, const Cube< Bool > &yinFlags, Int method, Bool goodIsTrue=False, Bool extrapolate=False)
Interpolate in the middle dimension of 3D array yin whose x coordinates along this dimension are give...
static void interpolatey(Cube< Range > &yout, const Vector< Domain > &xout, const Vector< Domain > &xin, const Cube< Range > &yin, Int method)
Interpolate in the middle axis in 3D array (yin) whose x coordinates along the this dimension are giv...
static void interpolatePtr(PtrBlock< Range * > &yout, PtrBlock< Bool * > &youtFlags, Int ny, const Vector< Domain > &xout, const Vector< Domain > &xin, const PtrBlock< const Range * > &yin, const PtrBlock< const Bool * > &yinFlags, Int method, Bool goodIsTrue, Bool extrapolate)
Interpolate the y-vectors of length ny from x values xin to xout.
static void interpolate(Array< Range > &yout, const Vector< Domain > &xout, const Vector< Domain > &xin, const Array< Range > &yin, Int method)
Interpolate in the last dimension of array yin whose x coordinates along this dimension are given by ...
static void interpolatePtr(PtrBlock< Range * > &yout, Int ny, const Vector< Domain > &xout, const Vector< Domain > &xin, const PtrBlock< const Range * > &yin, Int method)
Interpolate the y-vectors of length ny from x values xin to xout.
static void interpolate(Array< Range > &yout, const Block< Domain > &xout, const Block< Domain > &xin, const Array< Range > &yin, Int method)
deprecated version of previous function using Blocks - no longer needed now that Vector has a fast in...
static void interpolate(Array< Range > &yout, Array< Bool > &youtFlags, const Vector< Domain > &xout, const Vector< Domain > &xin, const Array< Range > &yin, const Array< Bool > &yinFlags, Int method, Bool goodIsTrue=False, Bool extrapolate=False)
Interpolate in the last dimension of array yin whose x coordinates along this dimension are given by ...
static void interpolate(Array< Range > &yout, Array< Bool > &youtFlags, const Block< Domain > &xout, const Block< Domain > &xin, const Array< Range > &yin, const Array< Bool > &yinFlags, Int method, Bool goodIsTrue=False, Bool extrapolate=False)
deprecated version of previous function using Blocks - no longer needed now that Vector has a fast in...
static void interpolateyPtr(PtrBlock< Range * > &yout, Int na, Int nb, Int nc, const Vector< Domain > &xout, const Vector< Domain > &xin, const PtrBlock< const Range * > &yin, Int method)
Interpolate along yaxis
static void polynomialInterpolation(PtrBlock< Range * > &yout, Int ny, const Vector< Domain > &xout, const Vector< Domain > &xin, const PtrBlock< const Range * > &yin, Int order)
Interpolate the y-vectors of length ny from x values xin to xout using polynomial interpolation with ...
A drop-in replacement for Block<T*>.
Definition: Block.h:814
this file contains all the compiler specific defines
Definition: mainpage.dox:28
const Bool False
Definition: aipstype.h:44
int Int
Definition: aipstype.h:50
bool Bool
Define the standard types used by Casacore.
Definition: aipstype.h:42