source: trunk/src/MathUtils.cpp@ 1645

Last change on this file since 1645 was 1570, checked in by Malte Marquarding, 15 years ago

Ticket #167: c++ part of running polynomial smoothing

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 7.2 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 <scimath/Fitting/LinearFit.h>
42#include <scimath/Functionals/Polynomial.h>
43#include <scimath/Mathematics/AutoDiff.h>
44
45
46#include "MathUtils.h"
47
48using namespace casa;
49
50float mathutil::statistics(const String& which,
51 const MaskedArray<Float>& data)
52{
53 String str(which);
54 str.upcase();
55 if (str.contains(String("MIN"))) {
56 return min(data);
57 } else if (str.contains(String("MAX"))) {
58 return max(data);
59 } else if (str.contains(String("SUMSQ"))) {
60 return sumsquares(data);
61 } else if (str.contains(String("SUM"))) {
62 return sum(data);
63 } else if (str.contains(String("MEAN"))) {
64 return mean(data);
65 } else if (str.contains(String("VAR"))) {
66 return variance(data);
67 } else if (str.contains(String("STDDEV"))) {
68 return stddev(data);
69 } else if (str.contains(String("AVDEV"))) {
70 return avdev(data);
71 } else if (str.contains(String("RMS"))) {
72 uInt n = data.nelementsValid();
73 return sqrt(sumsquares(data)/n);
74 } else if (str.contains(String("MED"))) {
75 return median(data);
76 }
77 return 0.0;
78}
79
80
81void mathutil::replaceMaskByZero(Vector<Float>& data, const Vector<Bool>& mask)
82{
83 for (uInt i=0; i<data.nelements(); i++) {
84 if (!mask[i]) data[i] = 0.0;
85 }
86}
87
88
89std::vector<std::string> mathutil::tovectorstring(const Vector<String>& in)
90{
91 std::vector<std::string> out;
92 out.reserve(in.nelements());
93 for (Array<String>::const_iterator it = in.begin(); it != in.end(); ++it) {
94 out.push_back(*it);
95 }
96 return out;
97}
98
99Vector<String> mathutil::toVectorString(const std::vector<std::string>& in)
100{
101 Vector<String> out(in.size());
102 Array<String>::iterator oit = out.begin();
103 for (std::vector<std::string>::const_iterator it=in.begin() ;
104 it != in.end(); ++it,++oit) {
105 *oit = *it;
106 }
107 return out;
108}
109
110void mathutil::hanning(Vector<Float>& out, Vector<Bool>& outmask,
111 const Vector<Float>& in, const Vector<Bool>& mask,
112 Bool relaxed, Bool ignoreOther) {
113 Vector< Vector<Float> > weights(8);
114 Vector<Float> vals(3);
115 vals = 0.0;weights[0] = vals;// FFF
116 vals[0] = 1.0; vals[1] = 0.0; vals[2] = 0.0; weights[1] = vals;// TFF
117 vals[0] = 0.0; vals[1] = 1.0; vals[2] = 0.0; weights[2] = vals;// FTF
118 vals[0] = 1.0/3.0; vals[1] = 2.0/3.0; vals[2] = 0.0; weights[3] = vals;// TTF
119 vals[0] = 0.0; vals[1] = 0.0; vals[2] = 1.0;weights[4] = vals;// FFT
120 vals[0] = 0.5; vals[1] = 0.0; vals[2] = 0.5; weights[5] = vals;// TFT
121 vals[0] = 0.0; vals[1] = 2.0/3.0; vals[2] = 1.0/3.0; weights[6] = vals;// FTT
122 vals[0] = 0.25; vals[1] = 0.5; vals[2] = 0.25; weights[7] = vals;// TTT
123 // Chris' case
124 Vector<Bool> weighted(8);
125 if (relaxed) {
126 weighted = False;
127 weighted[7] = True;
128
129 } else {
130 weighted = True;
131 weighted[0] = False;
132 }
133
134 out.resize(in.nelements());
135 outmask.resize(mask.nelements());
136 // make special case for first and last
137 /// ...here
138 // loop from 1..n-2
139 out.resize(in.nelements());
140 out[0] = in[0];out[out.nelements()-1] = in[in.nelements()-1];
141 outmask.resize(mask.nelements());
142 outmask = False;
143 uInt m;Vector<Float>* w;
144 for (uInt i=1; i < out.nelements()-1;++i) {
145 m = mask[i-1] + 2*mask[i] + 4*mask[i+1];
146 w = &(weights[m]);
147 if (weighted[m]) {
148 out[i] = (*w)[0]*in[i-1] + (*w)[1]*in[i] + (*w)[2]*in[i+1];
149 outmask[i] = True;
150 } else { // mask it
151 out[i] = in[i];//use arbitrary value
152 outmask[i] = False;
153 }
154 }
155}
156
157
158void mathutil::runningMedian(Vector<Float>& out, Vector<Bool>& outflag,
159 const Vector<Float>& in, const Vector<Bool>& flag,
160 float width)
161{
162 Int hwidth = Int(width+0.5);
163 Int fwidth = hwidth*2+1;
164 out.resize(in.nelements());
165 outflag.resize(flag.nelements());
166 MedianSlider ms(hwidth);
167 Slice sl(0, fwidth-1);
168 Float medval = ms.add(const_cast<Vector<Float>& >(in)(sl),
169 const_cast<Vector<Bool>& >(flag)(sl));
170 uInt n = in.nelements();
171 for (uInt i=hwidth; i<(n-hwidth); ++i) {
172 // add data value
173 out[i] = ms.add(in[i+hwidth], flag[i+hwidth]);
174 outflag[i] = (ms.nval() == 0);
175 }
176 // replicate edge values from first value with full width of values
177 for (uInt i=0;i<hwidth;++i) {
178 out[i] = out[hwidth];
179 outflag[i] = outflag[hwidth];
180 out[n-1-i] = out[n-1-hwidth];
181 outflag[n-1-i] = outflag[n-1-hwidth];
182 }
183}
184
185void mathutil::polyfit(Vector<Float>& out, Vector<Bool>& outmask,
186 const Vector<Float>& in, const Vector<Bool>& mask,
187 float width, int order)
188{
189 Int hwidth = Int(width+0.5);
190 Int fwidth = hwidth*2+1;
191 out.resize(in.nelements());
192 outmask.resize(mask.nelements());
193 LinearFit<Float> fitter;
194 Polynomial<Float> poly(order);
195 fitter.setFunction(poly);
196 Vector<Float> sigma(fwidth);
197 sigma = 1.0;
198 Vector<Float> parms;
199 Vector<Float> x(fwidth);
200 indgen(x);
201
202 uInt n = in.nelements();
203
204 for (uInt i=hwidth; i<(n-hwidth); ++i) {
205 // add data value
206 if (mask[i]) {
207 Slice sl(i-hwidth, fwidth);
208 const Vector<Float> &y = const_cast<Vector<Float>& >(in)(sl);
209 const Vector<Bool> &m = const_cast<Vector<Bool>& >(mask)(sl);
210 parms = fitter.fit(x, y, sigma, &m);
211
212 poly.setCoefficients(parms);
213 out[i] = poly(x[hwidth]);//cout << in[i] <<"->"<<out[i]<<endl;
214 } else {
215 out[i] = in[i];
216 }
217 outmask[i] = mask[i];
218 }
219 // replicate edge values from first value with full width of values
220 for (uInt i=0;i<hwidth;++i) {
221 out[i] = out[hwidth];
222 outmask[i] = outmask[hwidth];
223 out[n-1-i] = out[n-1-hwidth];
224 outmask[n-1-i] = outmask[n-1-hwidth];
225 }
226}
Note: See TracBrowser for help on using the repository browser.