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