source: trunk/src/MathUtils.cpp@ 1436

Last change on this file since 1436 was 1412, checked in by Malte Marquarding, 17 years ago

use stl type array iterators as in the new version of casacore

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 5.9 KB
Line 
1//#---------------------------------------------------------------------------
2//# MathUtilities.cc: General math operations
3//#---------------------------------------------------------------------------
4//# Copyright (C) 2004
5//# ATNF
6//#
7//# This program is free software; you can redistribute it and/or modify it
8//# under the terms of the GNU General Public License as published by the Free
9//# Software Foundation; either version 2 of the License, or (at your option)
10//# any later version.
11//#
12//# This program is distributed in the hope that it will be useful, but
13//# WITHOUT ANY WARRANTY; without even the implied warranty of
14//# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
15//# Public License for more details.
16//#
17//# You should have received a copy of the GNU General Public License along
18//# with this program; if not, write to the Free Software Foundation, Inc.,
19//# 675 Massachusetts Ave, Cambridge, MA 02139, USA.
20//#
21//# Correspondence concerning this software should be addressed as follows:
22//# Internet email: Malte.Marquarding@csiro.au
23//# Postal address: Malte Marquarding,
24//# Australia Telescope National Facility,
25//# P.O. Box 76,
26//# Epping, NSW, 2121,
27//# AUSTRALIA
28//#
29//# $Id:
30//#---------------------------------------------------------------------------
31
32#include <casa/aips.h>
33#include <casa/Arrays/Vector.h>
34#include <casa/Arrays/Slice.h>
35#include <casa/Arrays/MaskedArray.h>
36#include <casa/Arrays/MaskArrMath.h>
37#include <casa/Arrays/VectorSTLIterator.h>
38#include <casa/BasicSL/String.h>
39#include <scimath/Mathematics/MedianSlider.h>
40
41#include "MathUtils.h"
42
43using namespace casa;
44
45float mathutil::statistics(const String& which,
46 const MaskedArray<Float>& data)
47{
48 String str(which);
49 str.upcase();
50 if (str.contains(String("MIN"))) {
51 return min(data);
52 } else if (str.contains(String("MAX"))) {
53 return max(data);
54 } else if (str.contains(String("SUMSQ"))) {
55 return sumsquares(data);
56 } else if (str.contains(String("SUM"))) {
57 return sum(data);
58 } else if (str.contains(String("MEAN"))) {
59 return mean(data);
60 } else if (str.contains(String("VAR"))) {
61 return variance(data);
62 } else if (str.contains(String("STDDEV"))) {
63 return stddev(data);
64 } else if (str.contains(String("AVDEV"))) {
65 return avdev(data);
66 } else if (str.contains(String("RMS"))) {
67 uInt n = data.nelementsValid();
68 return sqrt(sumsquares(data)/n);
69 } else if (str.contains(String("MED"))) {
70 return median(data);
71 }
72 return 0.0;
73}
74
75
76void mathutil::replaceMaskByZero(Vector<Float>& data, const Vector<Bool>& mask)
77{
78 for (uInt i=0; i<data.nelements(); i++) {
79 if (!mask[i]) data[i] = 0.0;
80 }
81}
82
83
84std::vector<std::string> mathutil::tovectorstring(const Vector<String>& in)
85{
86 std::vector<std::string> out;
87 out.reserve(in.nelements());
88 for (Array<String>::const_iterator it = in.begin(); it != in.end(); ++it) {
89 out.push_back(*it);
90 }
91 return out;
92}
93
94Vector<String> mathutil::toVectorString(const std::vector<std::string>& in)
95{
96 Vector<String> out(in.size());
97 Array<String>::iterator oit = out.begin();
98 for (std::vector<std::string>::const_iterator it=in.begin() ;
99 it != in.end(); ++it,++oit) {
100 *oit = *it;
101 }
102 return out;
103}
104
105void mathutil::hanning(Vector<Float>& out, Vector<Bool>& outmask,
106 const Vector<Float>& in, const Vector<Bool>& mask,
107 Bool relaxed, Bool ignoreOther) {
108 Vector< Vector<Float> > weights(8);
109 Vector<Float> vals(3);
110 vals = 0.0;weights[0] = vals;// FFF
111 vals[0] = 1.0; vals[1] = 0.0; vals[2] = 0.0; weights[1] = vals;// TFF
112 vals[0] = 0.0; vals[1] = 1.0; vals[2] = 0.0; weights[2] = vals;// FTF
113 vals[0] = 1.0/3.0; vals[1] = 2.0/3.0; vals[2] = 0.0; weights[3] = vals;// TTF
114 vals[0] = 0.0; vals[1] = 0.0; vals[2] = 1.0;weights[4] = vals;// FFT
115 vals[0] = 0.5; vals[1] = 0.0; vals[2] = 0.5; weights[5] = vals;// TFT
116 vals[0] = 0.0; vals[1] = 2.0/3.0; vals[2] = 1.0/3.0; weights[6] = vals;// FTT
117 vals[0] = 0.25; vals[1] = 0.5; vals[2] = 0.25; weights[7] = vals;// TTT
118 // Chris' case
119 Vector<Bool> weighted(8);
120 if (relaxed) {
121 weighted = False;
122 weighted[7] = True;
123
124 } else {
125 weighted = True;
126 weighted[0] = False;
127 }
128
129 out.resize(in.nelements());
130 outmask.resize(mask.nelements());
131 // make special case for first and last
132 /// ...here
133 // loop from 1..n-2
134 out.resize(in.nelements());
135 out[0] = in[0];out[out.nelements()-1] = in[in.nelements()-1];
136 outmask.resize(mask.nelements());
137 outmask = False;
138 uInt m;Vector<Float>* w;
139 for (uInt i=1; i < out.nelements()-1;++i) {
140 m = mask[i-1] + 2*mask[i] + 4*mask[i+1];
141 w = &(weights[m]);
142 if (weighted[m]) {
143 out[i] = (*w)[0]*in[i-1] + (*w)[1]*in[i] + (*w)[2]*in[i+1];
144 outmask[i] = True;
145 } else { // mask it
146 out[i] = in[i];//use arbitrary value
147 outmask[i] = False;
148 }
149 }
150}
151
152
153void mathutil::runningMedian(Vector<Float>& out, Vector<Bool>& outflag,
154 const Vector<Float>& in, const Vector<Bool>& flag,
155 float width)
156{
157 Int hwidth = Int(width+0.5);
158 Int fwidth = hwidth*2+1;
159 out.resize(in.nelements());
160 outflag.resize(flag.nelements());
161 MedianSlider ms(hwidth);
162 Slice sl(0, fwidth-1);
163 Float medval = ms.add(const_cast<Vector<Float>& >(in)(sl),
164 const_cast<Vector<Bool>& >(flag)(sl));
165 uInt n = in.nelements();
166 for (uInt i=hwidth; i<(n-hwidth); ++i) {
167 // add data value
168 out[i] = ms.add(in[i+hwidth], flag[i+hwidth]);
169 outflag[i] = (ms.nval() == 0);
170 }
171 // replicate edge values from fisrt value with full width of values
172 for (uInt i=0;i<hwidth;++i) {
173 out[i] = out[hwidth];
174 outflag[i] = outflag[hwidth];
175 out[n-1-i] = out[n-1-hwidth];
176 outflag[n-1-i] = outflag[n-1-hwidth];
177 }
178}
Note: See TracBrowser for help on using the repository browser.