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$
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 
34 namespace casacore { //# NAMESPACE CASACORE - BEGIN
35 
36 // <summary>
37 // Mathematical and logical operations for Array parts.
38 // </summary>
39 // <reviewed reviewer="UNKNOWN" date="before2004/08/25" tests="tArray">
40 //
41 // <prerequisite>
42 // <li> <linkto class=Array>Array</linkto>
43 // </prerequisite>
44 //
45 // <etymology>
46 // This file contains global functions which perform part by part
47 // mathematical or logical operations on arrays.
48 // </etymology>
49 //
50 // <synopsis>
51 // These functions perform chunk by chunk mathematical operations on
52 // arrays.
53 // In particular boxed and sliding operations are possible. E.g. to calculate
54 // the median in sliding windows making it possible to subtract the background
55 // in an image.
56 //
57 // The operations to be performed are defined by means of functors that
58 // reduce an array subset to a scalar. Those functors are wrappers for
59 // ArrayMath and ArrayLogical functions like sum, median, and ntrue.
60 //
61 // The <src>partialXX</src> functions are a special case of the
62 // <src>BoxedArrayMath</src> function.
63 // They reduce one or more entire axes which can be done in a faster way than
64 // the more general <src>boxedArrayMath</src> function.
65 // </synopsis>
66 //
67 // <example>
68 // <srcblock>
69 // Array<Double> data(...);
70 // Array<Double> means = partialMeans (data, IPosition(2,0,1));
71 // </srcblock>
72 // This example calculates the mean of each plane in the data array.
73 // </example>
74 //
75 // <example>
76 // <srcblock>
77 // IPosition shp = data.shape();
78 // Array<Double> means = boxedArrayMath (data, IPosition(2,shp[0],shp[1]),
79 // SumFunc<Double>());
80 // </srcblock>
81 // does the same as the first example.
82 // Note that in this example the box is formed by the entire axes, but it
83 // could also be a subset of it to average, say, boxes of 5*5 elements.
84 // </example>
85 //
86 // <linkfrom anchor="Array mathematical operations" classes="Array Vector Matrix Cube">
87 // <here>Array mathematical operations</here> -- Mathematical operations for
88 // Arrays.
89 // </linkfrom>
90 //
91 // <group name="Array partial operations">
92 
93 
94 // Determine the sum, product, etc. for the given axes only.
95 // The result is an array with a shape formed by the remaining axes.
96 // For example, for an array with shape [3,4,5], collapsing axis 0
97 // results in an array with shape [4,5] containing, say, the sum for
98 // each X line.
99 // Summing for axes 0 and 2 results in an array with shape [4] containing,
100 // say, the sum for each XZ plane.
101 // <note>
102 // ArrayLogical.h contains the functions ntrue, nfalse, partialNTrue and
103 // partialNFalse to count the number of true or false elements in an array.
104 // </note>
105 // <group>
106 template<class T> Array<T> partialSums (const Array<T>& array,
107  const IPosition& collapseAxes);
108 template<class T> Array<T> partialProducts (const Array<T>& array,
109  const IPosition& collapseAxes);
110 template<class T> Array<T> partialMins (const Array<T>& array,
111  const IPosition& collapseAxes);
112 template<class T> Array<T> partialMaxs (const Array<T>& array,
113  const IPosition& collapseAxes);
114 template<class T> Array<T> partialMeans (const Array<T>& array,
115  const IPosition& collapseAxes);
116 template<class T> inline Array<T> partialVariances (const Array<T>& array,
117  const IPosition& collapseAxes)
118 {
119  return partialVariances (array, collapseAxes,
120  partialMeans (array, collapseAxes));
121 }
122 template<class T> Array<T> partialVariances (const Array<T>& array,
123  const IPosition& collapseAxes,
124  const Array<T>& means);
125 template<class T> inline Array<T> partialStddevs (const Array<T>& array,
126  const IPosition& collapseAxes)
127 {
128  return sqrt (partialVariances (array, collapseAxes,
129  partialMeans (array, collapseAxes)));
130 }
131 template<class T> inline Array<T> partialStddevs (const Array<T>& array,
132  const IPosition& collapseAxes,
133  const Array<T>& means)
134 {
135  return sqrt (partialVariances (array, collapseAxes, means));
136 }
137 template<class T> inline Array<T> partialAvdevs (const Array<T>& array,
138  const IPosition& collapseAxes)
139 {
140  return partialAvdevs (array, collapseAxes,
141  partialMeans (array, collapseAxes));
142 }
143 template<class T> Array<T> partialAvdevs (const Array<T>& array,
144  const IPosition& collapseAxes,
145  const Array<T>& means);
146 template<class T> Array<T> partialRmss (const Array<T>& array,
147  const IPosition& collapseAxes);
148 template<class T> Array<T> partialMedians (const Array<T>& array,
149  const IPosition& collapseAxes,
150  Bool takeEvenMean=False,
151  Bool inPlace=False);
152 template<class T> Array<T> partialMadfms (const Array<T>& array,
153  const IPosition& collapseAxes,
154  Bool takeEvenMean=False,
155  Bool inPlace=False);
156 template<class T> Array<T> partialFractiles (const Array<T>& array,
157  const IPosition& collapseAxes,
158  Float fraction,
159  Bool inPlace=False);
160 template<class T> Array<T> partialInterFractileRanges (const Array<T>& array,
161  const IPosition& collapseAxes,
162  Float fraction,
163  Bool inPlace=False);
164 template<class T> Array<T> partialInterHexileRanges (const Array<T>& array,
165  const IPosition& collapseAxes,
166  Bool inPlace=False)
167  { return partialInterFractileRanges (array, collapseAxes, 1./6., inPlace); }
168 template<class T> Array<T> partialInterQuartileRanges (const Array<T>& array,
169  const IPosition& collapseAxes,
170  Bool inPlace=False)
171  { return partialInterFractileRanges (array, collapseAxes, 0.25, inPlace); }
172 // </group>
173 
174 
175 
176 template<typename T> class SumFunc {
177 public:
178  T operator() (const Array<T>& arr) const { return sum(arr); }
179 };
180 template<typename T> class ProductFunc {
181 public:
182  T operator() (const Array<T>& arr) const { return product(arr); }
183 };
184 template<typename T> class MinFunc {
185 public:
186  T operator() (const Array<T>& arr) const { return min(arr); }
187 };
188 template<typename T> class MaxFunc {
189 public:
190  T operator() (const Array<T>& arr) const { return max(arr); }
191 };
192 template<typename T> class MeanFunc {
193 public:
194  T operator() (const Array<T>& arr) const { return mean(arr); }
195 };
196 template<typename T> class VarianceFunc {
197 public:
198  T operator() (const Array<T>& arr) const { return variance(arr); }
199 };
200 template<typename T> class StddevFunc {
201 public:
202  T operator() (const Array<T>& arr) const { return stddev(arr); }
203 };
204 template<typename T> class AvdevFunc {
205 public:
206  T operator() (const Array<T>& arr) const { return avdev(arr); }
207 };
208 template<typename T> class RmsFunc {
209 public:
210  T operator() (const Array<T>& arr) const { return rms(arr); }
211 };
212 template<typename T> class MedianFunc {
213 public:
214  explicit MedianFunc (Bool sorted=False, Bool takeEvenMean=True,
215  Bool inPlace = False)
216  : itsSorted(sorted), itsTakeEvenMean(takeEvenMean), itsInPlace(inPlace) {}
217  T operator() (const Array<T>& arr) const
218  { return median(arr, itsTmp, itsSorted, itsTakeEvenMean, itsInPlace); }
219 private:
223  mutable Block<T> itsTmp;
224 };
225 template<typename T> class MadfmFunc {
226 public:
227  explicit MadfmFunc(Bool sorted = False, Bool takeEvenMean = True,
228  Bool inPlace = False)
229  : itsSorted(sorted), itsTakeEvenMean(takeEvenMean), itsInPlace(inPlace) {}
230  T operator()(const Array<T>& arr) const
231  { return madfm(arr, itsTmp, itsSorted, itsTakeEvenMean, itsInPlace); }
232 private:
237 };
238 template<typename T> class FractileFunc {
239 public:
240  explicit FractileFunc (Float fraction,
241  Bool sorted = False, Bool inPlace = False)
242  : itsFraction(fraction), itsSorted(sorted), itsInPlace(inPlace) {}
243  T operator() (const Array<T>& arr) const
244  { return fractile(arr, itsTmp, itsFraction, itsSorted, itsInPlace); }
245 private:
246  float itsFraction;
249  mutable Block<T> itsTmp;
250 };
251 template<typename T> class InterFractileRangeFunc {
252 public:
253  explicit InterFractileRangeFunc(Float fraction,
254  Bool sorted = False, Bool inPlace = False)
255  : itsFraction(fraction), itsSorted(sorted), itsInPlace(inPlace) {}
256  T operator()(const Array<T>& arr) const
257  { return interFractileRange(arr, itsTmp, itsFraction,
258  itsSorted, itsInPlace); }
259 private:
260  float itsFraction;
264 };
265 template<typename T> class InterHexileRangeFunc: public InterFractileRangeFunc<T> {
266 public:
267  explicit InterHexileRangeFunc(Bool sorted = False, Bool inPlace = False)
268  : InterFractileRangeFunc<T> (1./6., sorted, inPlace)
269  {}
270 };
271 template<typename T> class InterQuartileRangeFunc: public InterFractileRangeFunc<T> {
272 public:
273  explicit InterQuartileRangeFunc(Bool sorted = False, Bool inPlace = False)
274  : InterFractileRangeFunc<T> (0.25, sorted, inPlace)
275  {}
276 };
277 
278 // Apply the given ArrayMath reduction function objects
279 // to each box in the array.
280 // <example>
281 // Downsample an array by taking the median of every [25,25] elements.
282 // <srcblock>
283 // Array<Float> downArr = boxedArrayMath(in, IPosition(2,25,25),
284 // MedianFunc<Float>());
285 // </srcblock>
286 // </example>
287 // The dimensionality of the array can be larger than the box; in that
288 // case the missing axes of the box are assumed to have length 1.
289 // A box axis length <= 0 means the full array axis.
290 template <typename T, typename FuncType>
291 Array<T> boxedArrayMath (const Array<T>& array,
292  const IPosition& boxSize,
293  const FuncType& funcObj);
294 
295 // Apply for each element in the array the given ArrayMath reduction function
296 // object to the box around that element. The full box is 2*halfBoxSize + 1.
297 // It can be used for arrays and boxes of any dimensionality; missing
298 // halfBoxSize values are set to 0.
299 // <example>
300 // Determine for each element in the array the median of a box
301 // with size [51,51] around that element:
302 // <srcblock>
303 // Array<Float> medians = slidingArrayMath(in, IPosition(2,25,25),
304 // MedianFunc<Float>());
305 // </srcblock>
306 // This is a potentially expensive operation. On a high-end PC it took
307 // appr. 27 seconds to get the medians for an array of [1000,1000] using
308 // a halfBoxSize of [50,50].
309 // </example>
310 // <br>The fillEdge argument determines how the edge is filled where
311 // no full boxes can be made. True means it is set to zero; False means
312 // that the edge is removed, thus the output array is smaller than the
313 // input array.
314 // <note> This brute-force method of determining the medians outperforms
315 // all kinds of smart implementations. For a vector it is about as fast
316 // as class <linkto class=MedianSlider>MedianSlider</linkto>, for a 2D array
317 // it is much, much faster.
318 // </note>
319 //# Do not use this template definition, as it is stricter. It precludes use
320 //# of double accumulators on float data.
321 //# template <typename T, template <typename T> class FuncType>
322 //# Array<T> slidingArrayMath (const Array<T>& array,
323 //# const IPosition& halfBoxSize,
324 //# const FuncType<T>& funcObj,
325 //# Bool fillEdge=True);
326 template <typename T, typename FuncType>
327 Array<T> slidingArrayMath (const Array<T>& array,
328  const IPosition& halfBoxSize,
329  const FuncType& funcObj,
330  Bool fillEdge=True);
331 
332 // </group>
333 
334 } //# NAMESPACE CASACORE - END
335 
336 #ifndef CASACORE_NO_AUTO_TEMPLATES
337 #include <casacore/casa/Arrays/ArrayPartMath.tcc>
338 #endif //# CASACORE_NO_AUTO_TEMPLATES
339 #endif
A Vector of integers, for indexing into Array<T> objects.
Definition: IPosition.h:119
TableExprNode means(const TableExprNode &array, const TableExprNodeSet &collapseAxes)
Definition: ExprNode.h:1799
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:1986
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)
MedianFunc(Bool sorted=False, Bool takeEvenMean=True, Bool inPlace=False)
LatticeExprNode min(const LatticeExprNode &left, const LatticeExprNode &right)
LatticeExprNode avdev(const LatticeExprNode &expr)
LatticeExprNode sqrt(const LatticeExprNode &expr)
bool Bool
Define the standard types used by Casacore.
Definition: aipstype.h:39
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:51
MadfmFunc(Bool sorted=False, Bool takeEvenMean=True, Bool inPlace=False)
const Bool False
Definition: aipstype.h:41
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)
LatticeExprNode mean(const LatticeExprNode &expr)
TableExprNode rms(const TableExprNode &array)
Definition: ExprNode.h:1737
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:40
this file contains all the compiler specific defines
Definition: mainpage.dox:28