ArrayPartMath.h

Go to the documentation of this file.
00001 //# ArrayPartMath.h: mathematics done on an array parts.
00002 //# Copyright (C) 1993,1994,1995,1996,1998,1999,2001,2003
00003 //# Associated Universities, Inc. Washington DC, USA.
00004 //# 
00005 //# This library is free software; you can redistribute it and/or modify it
00006 //# under the terms of the GNU Library General Public License as published by
00007 //# the Free Software Foundation; either version 2 of the License, or (at your
00008 //# option) any later version.
00009 //# 
00010 //# This library is distributed in the hope that it will be useful, but WITHOUT
00011 //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00012 //# FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Library General Public
00013 //# License for more details.
00014 //# 
00015 //# You should have received a copy of the GNU Library General Public License
00016 //# along with this library; if not, write to the Free Software Foundation,
00017 //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
00018 //# 
00019 //# Correspondence concerning AIPS++ should be addressed as follows:
00020 //#        Internet email: aips2-request@nrao.edu.
00021 //#        Postal address: AIPS++ Project Office
00022 //#                        National Radio Astronomy Observatory
00023 //#                        520 Edgemont Road
00024 //#                        Charlottesville, VA 22903-2475 USA
00025 //#
00026 //# $Id: ArrayPartMath.h 21262 2012-09-07 12:38:36Z gervandiepen $
00027 
00028 #ifndef CASA_ARRAYPARTMATH_H
00029 #define CASA_ARRAYPARTMATH_H
00030 
00031 #include <casacore/casa/aips.h>
00032 #include <casacore/casa/Arrays/ArrayMath.h>
00033 #include <casacore/casa/Arrays/ArrayMathBase.h>
00034 
00035 namespace casacore { //# NAMESPACE CASACORE - BEGIN
00036 
00037 // <summary>
00038 //    Mathematical and logical operations for Array parts.
00039 // </summary>
00040 // <reviewed reviewer="UNKNOWN" date="before2004/08/25" tests="tArray">
00041 //
00042 // <prerequisite>
00043 //   <li> <linkto class=Array>Array</linkto>
00044 // </prerequisite>
00045 //
00046 // <etymology>
00047 // This file contains global functions which perform part by part
00048 // mathematical or logical operations on arrays.
00049 // </etymology>
00050 //
00051 // <synopsis>
00052 // These functions perform chunk by chunk mathematical operations on
00053 // arrays.
00054 // In particular boxed and sliding operations are possible. E.g. to calculate
00055 // the median in sliding windows making it possible to subtract the background
00056 // in an image.
00057 //
00058 // The operations to be performed are defined by means of functors that
00059 // reduce an array subset to a scalar. Those functors are wrappers for
00060 // ArrayMath and ArrayLogical functions like sum, median, and ntrue. 
00061 //
00062 // The <src>partialXX</src> functions are a special case of the 
00063 // <src>BoxedArrayMath</src> function.
00064 // They reduce one or more entire axes which can be done in a faster way than
00065 // the more general <src>boxedArrayMath</src> function.
00066 // </synopsis>
00067 //
00068 // <example>
00069 // <srcblock>
00070 // Array<Double> data(...);
00071 // Array<Double> means = partialMeans (data, IPosition(2,0,1));
00072 // </srcblock>
00073 // This example calculates the mean of each plane in the data array.
00074 // </example>
00075 //
00076 // <example>
00077 // <srcblock>
00078 // IPosition shp = data.shape();
00079 // Array<Double> means = boxedArrayMath (data, IPosition(2,shp[0],shp[1]),
00080 //                                       SumFunc<Double>());
00081 // </srcblock>
00082 // does the same as the first example.
00083 // Note that in this example the box is formed by the entire axes, but it
00084 // could also be a subset of it to average, say, boxes of 5*5 elements.
00085 // </example>
00086 //
00087 // <linkfrom anchor="Array mathematical operations" classes="Array Vector Matrix Cube">
00088 //    <here>Array mathematical operations</here> -- Mathematical operations for
00089 //    Arrays.
00090 // </linkfrom>
00091 //
00092 // <group name="Array partial operations">
00093 
00094 
00095 // Determine the sum, product, etc. for the given axes only.
00096 // The result is an array with a shape formed by the remaining axes.
00097 // For example, for an array with shape [3,4,5], collapsing axis 0
00098 // results in an array with shape [4,5] containing, say, the sum for
00099 // each X line.
00100 // Summing for axes 0 and 2 results in an array with shape [4] containing,
00101 // say, the sum for each XZ plane.
00102 // <note>
00103 // ArrayLogical.h contains the functions ntrue, nfalse, partialNTrue and
00104 // partialNFalse to count the number of true or false elements in an array.
00105 // </note>
00106 // <group>
00107 template<class T> Array<T> partialSums (const Array<T>& array,
00108                                         const IPosition& collapseAxes);
00109 template<class T> Array<T> partialSumSqrs (const Array<T>& array,
00110                                            const IPosition& collapseAxes);
00111 template<class T> Array<T> partialProducts (const Array<T>& array,
00112                                             const IPosition& collapseAxes);
00113 template<class T> Array<T> partialMins (const Array<T>& array,
00114                                         const IPosition& collapseAxes);
00115 template<class T> Array<T> partialMaxs (const Array<T>& array,
00116                                         const IPosition& collapseAxes);
00117 template<class T> Array<T> partialMeans (const Array<T>& array,
00118                                          const IPosition& collapseAxes);
00119 template<class T> inline Array<T> partialVariances (const Array<T>& array,
00120                                              const IPosition& collapseAxes)
00121 {
00122     return partialVariances (array, collapseAxes,
00123                              partialMeans (array, collapseAxes));
00124 }
00125 template<class T> Array<T> partialVariances (const Array<T>& array,
00126                                              const IPosition& collapseAxes,
00127                                              const Array<T>& means);
00128 template<class T> inline Array<T> partialStddevs (const Array<T>& array,
00129                                            const IPosition& collapseAxes)
00130 {
00131     return sqrt (partialVariances (array, collapseAxes,
00132                                    partialMeans (array, collapseAxes)));
00133 }
00134 template<class T> inline Array<T> partialStddevs (const Array<T>& array,
00135                                            const IPosition& collapseAxes,
00136                                            const Array<T>& means)
00137 {
00138     return sqrt (partialVariances (array, collapseAxes, means));
00139 }
00140 template<class T> inline Array<T> partialAvdevs (const Array<T>& array,
00141                                           const IPosition& collapseAxes)
00142 {
00143     return partialAvdevs (array, collapseAxes,
00144                           partialMeans (array, collapseAxes));
00145 }
00146 template<class T> Array<T> partialAvdevs (const Array<T>& array,
00147                                           const IPosition& collapseAxes,
00148                                           const Array<T>& means);
00149 template<class T> Array<T> partialRmss (const Array<T>& array,
00150                                         const IPosition& collapseAxes);
00151 template<class T> Array<T> partialMedians (const Array<T>& array,
00152                                            const IPosition& collapseAxes,
00153                                            Bool takeEvenMean=False,
00154                                            Bool inPlace=False);
00155 template<class T> Array<T> partialMadfms (const Array<T>& array,
00156                                           const IPosition& collapseAxes,
00157                                           Bool takeEvenMean=False,
00158                                           Bool inPlace=False);
00159 template<class T> Array<T> partialFractiles (const Array<T>& array,
00160                                              const IPosition& collapseAxes,
00161                                              Float fraction,
00162                                              Bool inPlace=False);
00163 template<class T> Array<T> partialInterFractileRanges (const Array<T>& array,
00164                                                        const IPosition& collapseAxes,
00165                                                        Float fraction,
00166                                                        Bool inPlace=False);
00167 template<class T> Array<T> partialInterHexileRanges (const Array<T>& array,
00168                                                      const IPosition& collapseAxes,
00169                                                      Bool inPlace=False)
00170   { return partialInterFractileRanges (array, collapseAxes, 1./6., inPlace); }
00171 template<class T> Array<T> partialInterQuartileRanges (const Array<T>& array,
00172                                                       const IPosition& collapseAxes,
00173                                                       Bool inPlace=False)
00174   { return partialInterFractileRanges (array, collapseAxes, 0.25, inPlace); }
00175 // </group>
00176 
00177 
00178 
00179   // Define functors to perform a reduction function on an Array object.
00180   // Use virtual functions instead of templates to avoid code bloat
00181   // in partialArrayMath, etc.
00182   template<typename T> class SumFunc : public ArrayFunctorBase<T> {
00183   public:
00184     virtual ~SumFunc() {}
00185     virtual T operator() (const Array<T>& arr) const { return sum(arr); }
00186   };
00187   template<typename T> class SumSqrFunc : public ArrayFunctorBase<T> {
00188   public:
00189     virtual ~SumSqrFunc() {}
00190     virtual T operator() (const Array<T>& arr) const { return sumsqr(arr); }
00191   };
00192   template<typename T> class ProductFunc : public ArrayFunctorBase<T> {
00193   public:
00194     virtual ~ProductFunc() {}
00195     virtual T operator() (const Array<T>& arr) const { return product(arr); }
00196   };
00197   template<typename T> class MinFunc : public ArrayFunctorBase<T> {
00198   public:
00199     virtual ~MinFunc() {}
00200     virtual T operator() (const Array<T>& arr) const { return min(arr); }
00201   };
00202   template<typename T> class MaxFunc : public ArrayFunctorBase<T> {
00203   public:
00204     virtual ~MaxFunc() {}
00205     virtual T operator() (const Array<T>& arr) const { return max(arr); }
00206   };
00207   template<typename T> class MeanFunc : public ArrayFunctorBase<T> {
00208   public:
00209     virtual ~MeanFunc() {}
00210     virtual T operator() (const Array<T>& arr) const { return mean(arr); }
00211   };
00212   template<typename T> class VarianceFunc : public ArrayFunctorBase<T> {
00213   public:
00214     virtual ~VarianceFunc() {}
00215     virtual T operator() (const Array<T>& arr) const { return variance(arr); }
00216   };
00217   template<typename T> class StddevFunc : public ArrayFunctorBase<T> {
00218   public:
00219     virtual ~StddevFunc() {}
00220     virtual T operator() (const Array<T>& arr) const { return stddev(arr); }
00221   };
00222   template<typename T> class AvdevFunc : public ArrayFunctorBase<T> {
00223   public:
00224     virtual ~AvdevFunc() {}
00225     virtual T operator() (const Array<T>& arr) const { return avdev(arr); }
00226   };
00227   template<typename T> class RmsFunc : public ArrayFunctorBase<T> {
00228   public:
00229     virtual ~RmsFunc() {}
00230     virtual T operator() (const Array<T>& arr) const { return rms(arr); }
00231   };
00232   template<typename T> class MedianFunc : public ArrayFunctorBase<T> {
00233   public:
00234     explicit MedianFunc (Bool sorted=False, Bool takeEvenMean=True,
00235                           Bool inPlace = False)
00236       : itsSorted(sorted), itsTakeEvenMean(takeEvenMean), itsInPlace(inPlace) {}
00237     virtual ~MedianFunc() {}
00238     virtual T operator() (const Array<T>& arr) const
00239       { return median(arr, itsTmp, itsSorted, itsTakeEvenMean, itsInPlace); }
00240   private:
00241     Bool     itsSorted;
00242     Bool     itsTakeEvenMean;
00243     Bool     itsInPlace;
00244     mutable Block<T> itsTmp;
00245   };
00246   template<typename T> class MadfmFunc {
00247   public:
00248     explicit MadfmFunc(Bool sorted = False, Bool takeEvenMean = True,
00249                        Bool inPlace = False)
00250       : itsSorted(sorted), itsTakeEvenMean(takeEvenMean), itsInPlace(inPlace) {}
00251     virtual ~MadfmFunc() {}
00252     virtual T operator()(const Array<T>& arr) const
00253       { return madfm(arr, itsTmp, itsSorted, itsTakeEvenMean, itsInPlace); }
00254   private:
00255     Bool     itsSorted;
00256     Bool     itsTakeEvenMean;
00257     Bool     itsInPlace;
00258     mutable Block<Float> itsTmp;
00259   };
00260   template<typename T> class FractileFunc : public ArrayFunctorBase<T> {
00261   public:
00262     explicit FractileFunc (Float fraction,
00263                             Bool sorted = False, Bool inPlace = False)
00264       : itsFraction(fraction), itsSorted(sorted), itsInPlace(inPlace) {}
00265     virtual ~FractileFunc() {}
00266     virtual T operator() (const Array<T>& arr) const
00267       { return fractile(arr, itsTmp, itsFraction, itsSorted, itsInPlace); }
00268   private:
00269     float    itsFraction;
00270     Bool     itsSorted;
00271     Bool     itsInPlace;
00272     mutable Block<T> itsTmp;
00273   };
00274   template<typename T> class InterFractileRangeFunc {
00275   public:
00276     explicit InterFractileRangeFunc(Float fraction,
00277                                     Bool sorted = False, Bool inPlace = False)
00278       : itsFraction(fraction), itsSorted(sorted), itsInPlace(inPlace) {}
00279     virtual ~InterFractileRangeFunc() {}
00280     virtual T operator()(const Array<T>& arr) const
00281       { return interFractileRange(arr, itsTmp, itsFraction,
00282                                   itsSorted, itsInPlace); }
00283   private:
00284     float    itsFraction;
00285     Bool     itsSorted;
00286     Bool     itsInPlace;
00287     mutable Block<Float> itsTmp;
00288   };
00289   template<typename T> class InterHexileRangeFunc: public InterFractileRangeFunc<T> {
00290   public:
00291     explicit InterHexileRangeFunc(Bool sorted = False, Bool inPlace = False)
00292       : InterFractileRangeFunc<T> (1./6., sorted, inPlace)
00293     {}
00294     virtual ~InterHexileRangeFunc() {}
00295   };
00296   template<typename T> class InterQuartileRangeFunc: public InterFractileRangeFunc<T> {
00297   public:
00298     explicit InterQuartileRangeFunc(Bool sorted = False, Bool inPlace = False)
00299       : InterFractileRangeFunc<T> (0.25, sorted, inPlace)
00300     {} 
00301     virtual ~InterQuartileRangeFunc() {}
00302   };
00303 
00304 
00305 
00306   // Do partial reduction of an Array object. I.e., perform the operation
00307   // on a subset of the array axes (the collapse axes).
00308   template<typename T>
00309   inline Array<T> partialArrayMath (const Array<T>& a,
00310                                     const IPosition& collapseAxes,
00311                                     const ArrayFunctorBase<T>& funcObj)
00312   {
00313     Array<T> res;
00314     partialArrayMath (res, a, collapseAxes, funcObj);
00315     return res;
00316   }
00317   template<typename T, typename RES>
00318   void partialArrayMath (Array<RES>& res,
00319                          const Array<T>& a,
00320                          const IPosition& collapseAxes,
00321                          const ArrayFunctorBase<T,RES>& funcObj);
00322 
00323 
00324 // Apply the given ArrayMath reduction function objects
00325 // to each box in the array.
00326 // <example>
00327 // Downsample an array by taking the median of every [25,25] elements.
00328 // <srcblock>
00329 //    Array<Float> downArr = boxedArrayMath(in, IPosition(2,25,25),
00330 //                                          MedianFunc<Float>());
00331 // </srcblock>
00332 // </example>
00333 // The dimensionality of the array can be larger than the box; in that
00334 // case the missing axes of the box are assumed to have length 1.
00335 // A box axis length <= 0 means the full array axis.
00336   template<typename T>
00337   inline Array<T> boxedArrayMath (const Array<T>& a,
00338                                   const IPosition& boxSize,
00339                                   const ArrayFunctorBase<T>& funcObj)
00340   {
00341     Array<T> res;
00342     boxedArrayMath (res, a, boxSize, funcObj);
00343     return res;
00344   }
00345   template<typename T, typename RES>
00346   void boxedArrayMath (Array<RES>&,
00347                        const Array<T>& array,
00348                        const IPosition& boxSize,
00349                        const ArrayFunctorBase<T,RES>& funcObj);
00350 
00351 // Apply for each element in the array the given ArrayMath reduction function
00352 // object to the box around that element. The full box is 2*halfBoxSize + 1.
00353 // It can be used for arrays and boxes of any dimensionality; missing
00354 // halfBoxSize values are set to 0.
00355 // <example>
00356 // Determine for each element in the array the median of a box
00357 // with size [51,51] around that element:
00358 // <srcblock>
00359 //    Array<Float> medians = slidingArrayMath(in, IPosition(2,25,25),
00360 //                                            MedianFunc<Float>());
00361 // </srcblock>
00362 // This is a potentially expensive operation. On a high-end PC it took
00363 // appr. 27 seconds to get the medians for an array of [1000,1000] using
00364 // a halfBoxSize of [50,50].
00365 // </example>
00366 // <br>The fillEdge argument determines how the edge is filled where
00367 // no full boxes can be made. True means it is set to zero; False means
00368 // that the edge is removed, thus the output array is smaller than the
00369 // input array.
00370 // <note> This brute-force method of determining the medians outperforms
00371 // all kinds of smart implementations. For a vector it is about as fast
00372 // as class <linkto class=MedianSlider>MedianSlider</linkto>, for a 2D array
00373 // it is much, much faster.
00374 // </note>
00375   template<typename T>
00376   inline Array<T> slidingArrayMath (const Array<T>& a,
00377                                     const IPosition& halfBoxSize,
00378                                     const ArrayFunctorBase<T>& funcObj,
00379                                     Bool fillEdge=True)
00380   {
00381     Array<T> res;
00382     slidingArrayMath (res, a, halfBoxSize, funcObj, fillEdge);
00383     return res;
00384   }
00385   template<typename T, typename RES>
00386   void slidingArrayMath (Array<RES>& res,
00387                          const Array<T>& array,
00388                          const IPosition& halfBoxSize,
00389                          const ArrayFunctorBase<T,RES>& funcObj,
00390                          Bool fillEdge=True);
00391 
00392 // </group>
00393 
00394 // <group>
00395 // Helper functions for boxed and sliding functions.
00396 // Determine full box shape and shape of result for a boxed operation.
00397 void fillBoxedShape (const IPosition& shape, const IPosition& boxShape,
00398                      IPosition& fullBoxShape, IPosition& resultShape);
00399 // Determine the box end and shape of result for a sliding operation.
00400 // It returns False if the result is empty.
00401 Bool fillSlidingShape (const IPosition& shape, const IPosition& halfBoxSize,
00402                        IPosition& boxEnd, IPosition& resultShape);
00403 // </group>
00404 
00405 } //# NAMESPACE CASACORE - END
00406 
00407 #ifndef CASACORE_NO_AUTO_TEMPLATES
00408 #include <casacore/casa/Arrays/ArrayPartMath.tcc>
00409 #endif //# CASACORE_NO_AUTO_TEMPLATES
00410 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

Generated on 31 Aug 2016 for casa by  doxygen 1.6.1