casacore
ArrayPartMath.h
Go to the documentation of this file.
1 //# ArrayPartMath.h: mathematics done on an array parts.
2 //# Copyright (C) 1993,1994,1995,1996,1998,1999,2001,2003
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: ArrayPartMath.h 21262 2012-09-07 12:38:36Z gervandiepen $
27 
28 #ifndef CASA_ARRAYPARTMATH_H
29 #define CASA_ARRAYPARTMATH_H
30 
31 #include <casacore/casa/aips.h>
32 #include <casacore/casa/Arrays/ArrayMath.h>
33 #include <casacore/casa/Arrays/ArrayMathBase.h>
34 
35 namespace casacore { //# NAMESPACE CASACORE - BEGIN
36 
37 // <summary>
38 // Mathematical and logical operations for Array parts.
39 // </summary>
40 // <reviewed reviewer="UNKNOWN" date="before2004/08/25" tests="tArray">
41 //
42 // <prerequisite>
43 // <li> <linkto class=Array>Array</linkto>
44 // </prerequisite>
45 //
46 // <etymology>
47 // This file contains global functions which perform part by part
48 // mathematical or logical operations on arrays.
49 // </etymology>
50 //
51 // <synopsis>
52 // These functions perform chunk by chunk mathematical operations on
53 // arrays.
54 // In particular boxed and sliding operations are possible. E.g. to calculate
55 // the median in sliding windows making it possible to subtract the background
56 // in an image.
57 //
58 // The operations to be performed are defined by means of functors that
59 // reduce an array subset to a scalar. Those functors are wrappers for
60 // ArrayMath and ArrayLogical functions like sum, median, and ntrue.
61 //
62 // The <src>partialXX</src> functions are a special case of the
63 // <src>BoxedArrayMath</src> function.
64 // They reduce one or more entire axes which can be done in a faster way than
65 // the more general <src>boxedArrayMath</src> function.
66 // </synopsis>
67 //
68 // <example>
69 // <srcblock>
70 // Array<Double> data(...);
71 // Array<Double> means = partialMeans (data, IPosition(2,0,1));
72 // </srcblock>
73 // This example calculates the mean of each plane in the data array.
74 // </example>
75 //
76 // <example>
77 // <srcblock>
78 // IPosition shp = data.shape();
79 // Array<Double> means = boxedArrayMath (data, IPosition(2,shp[0],shp[1]),
80 // SumFunc<Double>());
81 // </srcblock>
82 // does the same as the first example.
83 // Note that in this example the box is formed by the entire axes, but it
84 // could also be a subset of it to average, say, boxes of 5*5 elements.
85 // </example>
86 //
87 // <linkfrom anchor="Array mathematical operations" classes="Array Vector Matrix Cube">
88 // <here>Array mathematical operations</here> -- Mathematical operations for
89 // Arrays.
90 // </linkfrom>
91 //
92 // <group name="Array partial operations">
93 
94 
95 // Determine the sum, product, etc. for the given axes only.
96 // The result is an array with a shape formed by the remaining axes.
97 // For example, for an array with shape [3,4,5], collapsing axis 0
98 // results in an array with shape [4,5] containing, say, the sum for
99 // each X line.
100 // Summing for axes 0 and 2 results in an array with shape [4] containing,
101 // say, the sum for each XZ plane.
102 // <note>
103 // ArrayLogical.h contains the functions ntrue, nfalse, partialNTrue and
104 // partialNFalse to count the number of true or false elements in an array.
105 // </note>
106 // <group>
107 template<class T> Array<T> partialSums (const Array<T>& array,
108  const IPosition& collapseAxes);
109 template<class T> Array<T> partialSumSqrs (const Array<T>& array,
110  const IPosition& collapseAxes);
111 template<class T> Array<T> partialProducts (const Array<T>& array,
112  const IPosition& collapseAxes);
113 template<class T> Array<T> partialMins (const Array<T>& array,
114  const IPosition& collapseAxes);
115 template<class T> Array<T> partialMaxs (const Array<T>& array,
116  const IPosition& collapseAxes);
117 template<class T> Array<T> partialMeans (const Array<T>& array,
118  const IPosition& collapseAxes);
119 template<class T> inline Array<T> partialVariances (const Array<T>& array,
120  const IPosition& collapseAxes)
121 {
122  return partialVariances (array, collapseAxes,
123  partialMeans (array, collapseAxes));
124 }
125 template<class T> Array<T> partialVariances (const Array<T>& array,
126  const IPosition& collapseAxes,
127  const Array<T>& means);
128 template<class T> inline Array<T> partialStddevs (const Array<T>& array,
129  const IPosition& collapseAxes)
130 {
131  return sqrt (partialVariances (array, collapseAxes,
132  partialMeans (array, collapseAxes)));
133 }
134 template<class T> inline Array<T> partialStddevs (const Array<T>& array,
135  const IPosition& collapseAxes,
136  const Array<T>& means)
137 {
138  return sqrt (partialVariances (array, collapseAxes, means));
139 }
140 template<class T> inline Array<T> partialAvdevs (const Array<T>& array,
141  const IPosition& collapseAxes)
142 {
143  return partialAvdevs (array, collapseAxes,
144  partialMeans (array, collapseAxes));
145 }
146 template<class T> Array<T> partialAvdevs (const Array<T>& array,
147  const IPosition& collapseAxes,
148  const Array<T>& means);
149 template<class T> Array<T> partialRmss (const Array<T>& array,
150  const IPosition& collapseAxes);
151 template<class T> Array<T> partialMedians (const Array<T>& array,
152  const IPosition& collapseAxes,
153  Bool takeEvenMean=False,
154  Bool inPlace=False);
155 template<class T> Array<T> partialMadfms (const Array<T>& array,
156  const IPosition& collapseAxes,
157  Bool takeEvenMean=False,
158  Bool inPlace=False);
159 template<class T> Array<T> partialFractiles (const Array<T>& array,
160  const IPosition& collapseAxes,
161  Float fraction,
162  Bool inPlace=False);
163 template<class T> Array<T> partialInterFractileRanges (const Array<T>& array,
164  const IPosition& collapseAxes,
165  Float fraction,
166  Bool inPlace=False);
167 template<class T> Array<T> partialInterHexileRanges (const Array<T>& array,
168  const IPosition& collapseAxes,
169  Bool inPlace=False)
170  { return partialInterFractileRanges (array, collapseAxes, 1./6., inPlace); }
171 template<class T> Array<T> partialInterQuartileRanges (const Array<T>& array,
172  const IPosition& collapseAxes,
173  Bool inPlace=False)
174  { return partialInterFractileRanges (array, collapseAxes, 0.25, inPlace); }
175 // </group>
176 
177 
178 
179  // Define functors to perform a reduction function on an Array object.
180  // Use virtual functions instead of templates to avoid code bloat
181  // in partialArrayMath, etc.
182  template<typename T> class SumFunc : public ArrayFunctorBase<T> {
183  public:
184  virtual ~SumFunc() {}
185  virtual T operator() (const Array<T>& arr) const { return sum(arr); }
186  };
187  template<typename T> class SumSqrFunc : public ArrayFunctorBase<T> {
188  public:
189  virtual ~SumSqrFunc() {}
190  virtual T operator() (const Array<T>& arr) const { return sumsqr(arr); }
191  };
192  template<typename T> class ProductFunc : public ArrayFunctorBase<T> {
193  public:
194  virtual ~ProductFunc() {}
195  virtual T operator() (const Array<T>& arr) const { return product(arr); }
196  };
197  template<typename T> class MinFunc : public ArrayFunctorBase<T> {
198  public:
199  virtual ~MinFunc() {}
200  virtual T operator() (const Array<T>& arr) const { return min(arr); }
201  };
202  template<typename T> class MaxFunc : public ArrayFunctorBase<T> {
203  public:
204  virtual ~MaxFunc() {}
205  virtual T operator() (const Array<T>& arr) const { return max(arr); }
206  };
207  template<typename T> class MeanFunc : public ArrayFunctorBase<T> {
208  public:
209  virtual ~MeanFunc() {}
210  virtual T operator() (const Array<T>& arr) const { return mean(arr); }
211  };
212  template<typename T> class VarianceFunc : public ArrayFunctorBase<T> {
213  public:
214  virtual ~VarianceFunc() {}
215  virtual T operator() (const Array<T>& arr) const { return variance(arr); }
216  };
217  template<typename T> class StddevFunc : public ArrayFunctorBase<T> {
218  public:
219  virtual ~StddevFunc() {}
220  virtual T operator() (const Array<T>& arr) const { return stddev(arr); }
221  };
222  template<typename T> class AvdevFunc : public ArrayFunctorBase<T> {
223  public:
224  virtual ~AvdevFunc() {}
225  virtual T operator() (const Array<T>& arr) const { return avdev(arr); }
226  };
227  template<typename T> class RmsFunc : public ArrayFunctorBase<T> {
228  public:
229  virtual ~RmsFunc() {}
230  virtual T operator() (const Array<T>& arr) const { return rms(arr); }
231  };
232  template<typename T> class MedianFunc : public ArrayFunctorBase<T> {
233  public:
234  explicit MedianFunc (Bool sorted=False, Bool takeEvenMean=True,
235  Bool inPlace = False)
236  : itsSorted(sorted), itsTakeEvenMean(takeEvenMean), itsInPlace(inPlace) {}
237  virtual ~MedianFunc() {}
238  virtual T operator() (const Array<T>& arr) const
239  { return median(arr, itsTmp, itsSorted, itsTakeEvenMean, itsInPlace); }
240  private:
244  mutable Block<T> itsTmp;
245  };
246  template<typename T> class MadfmFunc {
247  public:
248  explicit MadfmFunc(Bool sorted = False, Bool takeEvenMean = True,
249  Bool inPlace = False)
250  : itsSorted(sorted), itsTakeEvenMean(takeEvenMean), itsInPlace(inPlace) {}
251  virtual ~MadfmFunc() {}
252  virtual T operator()(const Array<T>& arr) const
253  { return madfm(arr, itsTmp, itsSorted, itsTakeEvenMean, itsInPlace); }
254  private:
259  };
260  template<typename T> class FractileFunc : public ArrayFunctorBase<T> {
261  public:
262  explicit FractileFunc (Float fraction,
263  Bool sorted = False, Bool inPlace = False)
264  : itsFraction(fraction), itsSorted(sorted), itsInPlace(inPlace) {}
265  virtual ~FractileFunc() {}
266  virtual T operator() (const Array<T>& arr) const
267  { return fractile(arr, itsTmp, itsFraction, itsSorted, itsInPlace); }
268  private:
269  float itsFraction;
272  mutable Block<T> itsTmp;
273  };
274  template<typename T> class InterFractileRangeFunc {
275  public:
276  explicit InterFractileRangeFunc(Float fraction,
277  Bool sorted = False, Bool inPlace = False)
278  : itsFraction(fraction), itsSorted(sorted), itsInPlace(inPlace) {}
280  virtual T operator()(const Array<T>& arr) const
281  { return interFractileRange(arr, itsTmp, itsFraction,
282  itsSorted, itsInPlace); }
283  private:
284  float itsFraction;
288  };
289  template<typename T> class InterHexileRangeFunc: public InterFractileRangeFunc<T> {
290  public:
291  explicit InterHexileRangeFunc(Bool sorted = False, Bool inPlace = False)
292  : InterFractileRangeFunc<T> (1./6., sorted, inPlace)
293  {}
295  };
296  template<typename T> class InterQuartileRangeFunc: public InterFractileRangeFunc<T> {
297  public:
298  explicit InterQuartileRangeFunc(Bool sorted = False, Bool inPlace = False)
299  : InterFractileRangeFunc<T> (0.25, sorted, inPlace)
300  {}
302  };
303 
304 
305 
306  // Do partial reduction of an Array object. I.e., perform the operation
307  // on a subset of the array axes (the collapse axes).
308  template<typename T>
310  const IPosition& collapseAxes,
311  const ArrayFunctorBase<T>& funcObj)
312  {
313  Array<T> res;
314  partialArrayMath (res, a, collapseAxes, funcObj);
315  return res;
316  }
317  template<typename T, typename RES>
318  void partialArrayMath (Array<RES>& res,
319  const Array<T>& a,
320  const IPosition& collapseAxes,
321  const ArrayFunctorBase<T,RES>& funcObj);
322 
323 
324 // Apply the given ArrayMath reduction function objects
325 // to each box in the array.
326 // <example>
327 // Downsample an array by taking the median of every [25,25] elements.
328 // <srcblock>
329 // Array<Float> downArr = boxedArrayMath(in, IPosition(2,25,25),
330 // MedianFunc<Float>());
331 // </srcblock>
332 // </example>
333 // The dimensionality of the array can be larger than the box; in that
334 // case the missing axes of the box are assumed to have length 1.
335 // A box axis length <= 0 means the full array axis.
336  template<typename T>
337  inline Array<T> boxedArrayMath (const Array<T>& a,
338  const IPosition& boxSize,
339  const ArrayFunctorBase<T>& funcObj)
340  {
341  Array<T> res;
342  boxedArrayMath (res, a, boxSize, funcObj);
343  return res;
344  }
345  template<typename T, typename RES>
346  void boxedArrayMath (Array<RES>&,
347  const Array<T>& array,
348  const IPosition& boxSize,
349  const ArrayFunctorBase<T,RES>& funcObj);
350 
351 // Apply for each element in the array the given ArrayMath reduction function
352 // object to the box around that element. The full box is 2*halfBoxSize + 1.
353 // It can be used for arrays and boxes of any dimensionality; missing
354 // halfBoxSize values are set to 0.
355 // <example>
356 // Determine for each element in the array the median of a box
357 // with size [51,51] around that element:
358 // <srcblock>
359 // Array<Float> medians = slidingArrayMath(in, IPosition(2,25,25),
360 // MedianFunc<Float>());
361 // </srcblock>
362 // This is a potentially expensive operation. On a high-end PC it took
363 // appr. 27 seconds to get the medians for an array of [1000,1000] using
364 // a halfBoxSize of [50,50].
365 // </example>
366 // <br>The fillEdge argument determines how the edge is filled where
367 // no full boxes can be made. True means it is set to zero; False means
368 // that the edge is removed, thus the output array is smaller than the
369 // input array.
370 // <note> This brute-force method of determining the medians outperforms
371 // all kinds of smart implementations. For a vector it is about as fast
372 // as class <linkto class=MedianSlider>MedianSlider</linkto>, for a 2D array
373 // it is much, much faster.
374 // </note>
375  template<typename T>
377  const IPosition& halfBoxSize,
378  const ArrayFunctorBase<T>& funcObj,
379  Bool fillEdge=True)
380  {
381  Array<T> res;
382  slidingArrayMath (res, a, halfBoxSize, funcObj, fillEdge);
383  return res;
384  }
385  template<typename T, typename RES>
386  void slidingArrayMath (Array<RES>& res,
387  const Array<T>& array,
388  const IPosition& halfBoxSize,
389  const ArrayFunctorBase<T,RES>& funcObj,
390  Bool fillEdge=True);
391 
392 // </group>
393 
394 // <group>
395 // Helper functions for boxed and sliding functions.
396 // Determine full box shape and shape of result for a boxed operation.
397 void fillBoxedShape (const IPosition& shape, const IPosition& boxShape,
398  IPosition& fullBoxShape, IPosition& resultShape);
399 // Determine the box end and shape of result for a sliding operation.
400 // It returns False if the result is empty.
401 Bool fillSlidingShape (const IPosition& shape, const IPosition& halfBoxSize,
402  IPosition& boxEnd, IPosition& resultShape);
403 // </group>
404 
405 } //# NAMESPACE CASACORE - END
406 
407 #ifndef CASACORE_NO_AUTO_TEMPLATES
408 #include <casacore/casa/Arrays/ArrayPartMath.tcc>
409 #endif //# CASACORE_NO_AUTO_TEMPLATES
410 #endif
A Vector of integers, for indexing into Array<T> objects.
Definition: IPosition.h:119
Define functors to perform a reduction function on an Array object.
TableExprNode means(const TableExprNode &array, const TableExprNodeSet &collapseAxes)
Definition: ExprNode.h:1918
LatticeExprNode median(const LatticeExprNode &expr)
T product(const TableVector< T > &tv)
Definition: TabVecMath.h:385
TableExprNode array(const TableExprNode &values, const TableExprNodeSet &shape)
Create an array of the given shape and fill it with the values.
Definition: ExprNode.h:2105
LatticeExprNode sum(const LatticeExprNode &expr)
LatticeExprNode max(const LatticeExprNode &left, const LatticeExprNode &right)
Array< T > partialInterQuartileRanges(const Array< T > &array, const IPosition &collapseAxes, Bool inPlace=False)
LatticeExprNode fractile(const LatticeExprNode &expr, const LatticeExprNode &fraction)
Determine the value of the element at the part fraction from the beginning of the given lattice...
InterFractileRangeFunc(Float fraction, Bool sorted=False, Bool inPlace=False)
void fillBoxedShape(const IPosition &shape, const IPosition &boxShape, IPosition &fullBoxShape, IPosition &resultShape)
Helper functions for boxed and sliding functions.
MedianFunc(Bool sorted=False, Bool takeEvenMean=True, Bool inPlace=False)
LatticeExprNode min(const LatticeExprNode &left, const LatticeExprNode &right)
LatticeExprNode avdev(const LatticeExprNode &expr)
Array< T > boxedArrayMath(const Array< T > &a, const IPosition &boxSize, const ArrayFunctorBase< T > &funcObj)
Apply the given ArrayMath reduction function objects to each box in the array.
Array< T > partialArrayMath(const Array< T > &a, const IPosition &collapseAxes, const ArrayFunctorBase< T > &funcObj)
Do partial reduction of an Array object.
Bool fillSlidingShape(const IPosition &shape, const IPosition &halfBoxSize, IPosition &boxEnd, IPosition &resultShape)
Determine the box end and shape of result for a sliding operation.
LatticeExprNode sqrt(const LatticeExprNode &expr)
Basic class for math on Array objects.
Definition: ArrayMathBase.h:57
bool Bool
Define the standard types used by Casacore.
Definition: aipstype.h:42
MaskedArray< T > boxedArrayMath(const MaskedArray< T > &array, const IPosition &boxSize, const FuncType &funcObj)
Apply the given ArrayMath reduction function objects to each box in the array.
LatticeExprNode stddev(const LatticeExprNode &expr)
float Float
Definition: aipstype.h:54
MadfmFunc(Bool sorted=False, Bool takeEvenMean=True, Bool inPlace=False)
const Bool False
Definition: aipstype.h:44
TableExprNode shape(const TableExprNode &array)
Function operating on any scalar or array resulting in a Double array containing the shape...
Definition: ExprNode.h:2163
template <class T, class U> class vector;
Definition: Array.h:169
FractileFunc(Float fraction, Bool sorted=False, Bool inPlace=False)
simple 1-D array
Definition: ArrayIO.h:47
Array< T > partialAvdevs(const Array< T > &array, const IPosition &collapseAxes)
Array< T > slidingArrayMath(const Array< T > &a, const IPosition &halfBoxSize, const ArrayFunctorBase< T > &funcObj, Bool fillEdge=True)
Apply for each element in the array the given ArrayMath reduction function object to the box around t...
LatticeExprNode mean(const LatticeExprNode &expr)
TableExprNode rms(const TableExprNode &array)
Definition: ExprNode.h:1856
Array< T > partialVariances(const Array< T > &array, const IPosition &collapseAxes)
Array< T > slidingArrayMath(const MaskedArray< T > &array, const IPosition &halfBoxSize, const FuncType &funcObj, Bool fillEdge=True)
Apply for each element in the array the given ArrayMath reduction function object to the box around t...
Array< T > partialInterHexileRanges(const Array< T > &array, const IPosition &collapseAxes, Bool inPlace=False)
LatticeExprNode variance(const LatticeExprNode &expr)
Array< T > partialStddevs(const Array< T > &array, const IPosition &collapseAxes, const Array< T > &means)
Array< T > partialStddevs(const Array< T > &array, const IPosition &collapseAxes)
const Bool True
Definition: aipstype.h:43
this file contains all the compiler specific defines
Definition: mainpage.dox:28