source: branches/alma/src/MathUtils.cpp@ 1435

Last change on this file since 1435 was 1400, checked in by TakTsutsumi, 17 years ago

added unflag parameter to flag function

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 5.8 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/BasicSL/String.h>
38#include <scimath/Mathematics/MedianSlider.h>
39
40#include "MathUtils.h"
41
42using namespace casa;
43
44float mathutil::statistics(const String& which,
45 const MaskedArray<Float>& data)
46{
47 String str(which);
48 str.upcase();
49 if (str.contains(String("MIN"))) {
50 return min(data);
51 } else if (str.contains(String("MAX"))) {
52 return max(data);
53 } else if (str.contains(String("SUMSQ"))) {
54 return sumsquares(data);
55 } else if (str.contains(String("SUM"))) {
56 return sum(data);
57 } else if (str.contains(String("MEAN"))) {
58 return mean(data);
59 } else if (str.contains(String("VAR"))) {
60 return variance(data);
61 } else if (str.contains(String("STDDEV"))) {
62 return stddev(data);
63 } else if (str.contains(String("AVDEV"))) {
64 return avdev(data);
65 } else if (str.contains(String("RMS"))) {
66 uInt n = data.nelementsValid();
67 return sqrt(sumsquares(data)/n);
68 } else if (str.contains(String("MED"))) {
69 return median(data);
70 }
71 return 0.0;
72}
73
74
75void mathutil::replaceMaskByZero(Vector<Float>& data, const Vector<Bool>& mask)
76{
77 for (uInt i=0; i<data.nelements(); i++) {
78 if (!mask[i]) data[i] = 0.0;
79 }
80}
81
82
83std::vector<std::string> mathutil::tovectorstring(const Vector<String>& in)
84{
85 std::vector<std::string> out;
86 Vector<String>::const_iterator it = in.begin();
87 for (uInt i=0; it != in.end(); ++it,++i) {
88 out.push_back(*it);
89 }
90 return out;
91}
92
93Vector<String> mathutil::toVectorString(const std::vector<std::string>& in)
94{
95 Vector<String> out(in.size());
96 uInt i=0;
97 std::vector<std::string>::const_iterator it;
98 for (it=in.begin();it != in.end();++it,++i) {
99 out[i] = casa::String(*it);
100 }
101 return out;
102}
103
104void mathutil::hanning(Vector<Float>& out, Vector<Bool>& outmask,
105 const Vector<Float>& in, const Vector<Bool>& mask,
106 Bool relaxed, Bool ignoreOther) {
107 Vector< Vector<Float> > weights(8);
108 Vector<Float> vals(3);
109 vals = 0.0;weights[0] = vals;// FFF
110 vals[0] = 1.0; vals[1] = 0.0; vals[2] = 0.0; weights[1] = vals;// TFF
111 vals[0] = 0.0; vals[1] = 1.0; vals[2] = 0.0; weights[2] = vals;// FTF
112 vals[0] = 1.0/3.0; vals[1] = 2.0/3.0; vals[2] = 0.0; weights[3] = vals;// TTF
113 vals[0] = 0.0; vals[1] = 0.0; vals[2] = 1.0;weights[4] = vals;// FFT
114 vals[0] = 0.5; vals[1] = 0.0; vals[2] = 0.5; weights[5] = vals;// TFT
115 vals[0] = 0.0; vals[1] = 2.0/3.0; vals[2] = 1.0/3.0; weights[6] = vals;// FTT
116 vals[0] = 0.25; vals[1] = 0.5; vals[2] = 0.25; weights[7] = vals;// TTT
117 // Chris' case
118 Vector<Bool> weighted(8);
119 if (relaxed) {
120 weighted = False;
121 weighted[7] = True;
122
123 } else {
124 weighted = True;
125 weighted[0] = False;
126 }
127
128 out.resize(in.nelements());
129 outmask.resize(mask.nelements());
130 // make special case for first and last
131 /// ...here
132 // loop from 1..n-2
133 out.resize(in.nelements());
134 out[0] = in[0];out[out.nelements()-1] = in[in.nelements()-1];
135 outmask.resize(mask.nelements());
136 outmask = False;
137 uInt m;Vector<Float>* w;
138 for (uInt i=1; i < out.nelements()-1;++i) {
139 m = mask[i-1] + 2*mask[i] + 4*mask[i+1];
140 w = &(weights[m]);
141 if (weighted[m]) {
142 out[i] = (*w)[0]*in[i-1] + (*w)[1]*in[i] + (*w)[2]*in[i+1];
143 outmask[i] = True;
144 } else { // mask it
145 out[i] = in[i];//use arbitrary value
146 outmask[i] = False;
147 }
148 }
149}
150
151
152void mathutil::runningMedian(Vector<Float>& out, Vector<Bool>& outflag,
153 const Vector<Float>& in, const Vector<Bool>& flag,
154 float width)
155{
156 Int hwidth = Int(width+0.5);
157 Int fwidth = hwidth*2+1;
158 out.resize(in.nelements());
159 outflag.resize(flag.nelements());
160 MedianSlider ms(hwidth);
161 Slice sl(0, fwidth-1);
162 Float medval = ms.add(const_cast<Vector<Float>& >(in)(sl),
163 const_cast<Vector<Bool>& >(flag)(sl));
164 uInt n = in.nelements();
165 for (uInt i=hwidth; i<(n-hwidth); ++i) {
166 // add data value
167 out[i] = ms.add(in[i+hwidth], flag[i+hwidth]);
168 outflag[i] = (ms.nval() == 0);
169 }
170 // replicate edge values from fisrt value with full width of values
171 for (uInt i=0;i<hwidth;++i) {
172 out[i] = out[hwidth];
173 outflag[i] = outflag[hwidth];
174 out[n-1-i] = out[n-1-hwidth];
175 outflag[n-1-i] = outflag[n-1-hwidth];
176 }
177}
Note: See TracBrowser for help on using the repository browser.