source: trunk/src/MathUtils.cpp@ 2687

Last change on this file since 2687 was 2670, checked in by Malte Marquarding, 12 years ago

Fix for #270: rmedian flags were incorrectly handled. Also remove excessive conversion step

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 10.3 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#include <casa/Exceptions/Error.h>
41
42#include <scimath/Fitting/LinearFit.h>
43#include <scimath/Functionals/Polynomial.h>
44#include <scimath/Mathematics/AutoDiff.h>
45
46#include <ctime>
47#include <sys/time.h>
48
49#include "MathUtils.h"
50
51using namespace casa;
52
53float mathutil::statistics(const String& which,
54 const MaskedArray<Float>& data)
55{
56 String str(which);
57 str.upcase();
58 if (str.matches(String("MIN"))) {
59 return min(data);
60 } else if (str.matches(String("MAX"))) {
61 return max(data);
62 } else if (str.matches(String("SUMSQ"))) {
63 return sumsquares(data);
64 } else if (str.matches(String("SUM"))) {
65 return sum(data);
66 } else if (str.matches(String("MEAN"))) {
67 return mean(data);
68 } else if (str.matches(String("VAR"))) {
69 return variance(data);
70 } else if (str.matches(String("STDDEV"))) {
71 return stddev(data);
72 } else if (str.matches(String("AVDEV"))) {
73 return avdev(data);
74 } else if (str.matches(String("RMS"))) {
75 uInt n = data.nelementsValid();
76 return sqrt(sumsquares(data)/n);
77 } else if (str.matches(String("MEDIAN"))) {
78 return median(data);
79 } else {
80 String msg = str + " is not a valid type of statistics";
81 throw(AipsError(msg));
82 }
83 return 0.0;
84}
85
86IPosition mathutil::minMaxPos(const String& which,
87 const MaskedArray<Float>& data)
88{
89 Float minVal, maxVal;
90 IPosition minPos(data.ndim(), 0), maxPos(data.ndim(), 0);
91 minMax(minVal, maxVal, minPos, maxPos, data);
92 String str(which);
93 str.upcase();
94 if (str.contains(String("MIN"))) {
95 return minPos;
96 } else if (str.contains(String("MAX"))) {
97 return maxPos;
98 } else {
99 String msg = str + " is not a valid type of statistics";
100 throw(AipsError(msg));
101 }
102 //return 0.0;
103}
104
105void mathutil::replaceMaskByZero(Vector<Float>& data, const Vector<Bool>& mask)
106{
107 for (uInt i=0; i<data.nelements(); i++) {
108 if (!mask[i]) data[i] = 0.0;
109 }
110}
111
112
113std::vector<std::string> mathutil::tovectorstring(const Vector<String>& in)
114{
115 std::vector<std::string> out;
116 out.reserve(in.nelements());
117 for (Array<String>::const_iterator it = in.begin(); it != in.end(); ++it) {
118 out.push_back(*it);
119 }
120 return out;
121}
122
123Vector<String> mathutil::toVectorString(const std::vector<std::string>& in)
124{
125 Vector<String> out(in.size());
126 Array<String>::iterator oit = out.begin();
127 for (std::vector<std::string>::const_iterator it=in.begin() ;
128 it != in.end(); ++it,++oit) {
129 *oit = *it;
130 }
131 return out;
132}
133
134void mathutil::hanning(Vector<Float>& out, Vector<Bool>& outmask,
135 const Vector<Float>& in, const Vector<Bool>& mask,
136 Bool relaxed, Bool ignoreOther) {
137 (void) ignoreOther; //suppress unused warning
138 Vector< Vector<Float> > weights(8);
139 Vector<Float> vals(3);
140 vals = 0.0;weights[0] = vals;// FFF
141 vals[0] = 1.0; vals[1] = 0.0; vals[2] = 0.0; weights[1] = vals;// TFF
142 vals[0] = 0.0; vals[1] = 1.0; vals[2] = 0.0; weights[2] = vals;// FTF
143 vals[0] = 1.0/3.0; vals[1] = 2.0/3.0; vals[2] = 0.0; weights[3] = vals;// TTF
144 vals[0] = 0.0; vals[1] = 0.0; vals[2] = 1.0;weights[4] = vals;// FFT
145 vals[0] = 0.5; vals[1] = 0.0; vals[2] = 0.5; weights[5] = vals;// TFT
146 vals[0] = 0.0; vals[1] = 2.0/3.0; vals[2] = 1.0/3.0; weights[6] = vals;// FTT
147 vals[0] = 0.25; vals[1] = 0.5; vals[2] = 0.25; weights[7] = vals;// TTT
148 // Chris' case
149 Vector<Bool> weighted(8);
150 if (relaxed) {
151 weighted = False;
152 weighted[7] = True;
153
154 } else {
155 weighted = True;
156 weighted[0] = False;
157 }
158
159 out.resize(in.nelements());
160 outmask.resize(mask.nelements());
161 // make special case for first and last
162 /// ...here
163 // loop from 1..n-2
164 out.resize(in.nelements());
165 out[0] = in[0];out[out.nelements()-1] = in[in.nelements()-1];
166 outmask.resize(mask.nelements());
167 outmask[0] = mask[0]; outmask[outmask.nelements()-1] = mask[mask.nelements()-1];
168 uInt m;Vector<Float>* w;
169 for (uInt i=1; i < out.nelements()-1;++i) {
170 m = mask[i-1] + 2*mask[i] + 4*mask[i+1];
171 w = &(weights[m]);
172 if (weighted[m]) {
173 out[i] = (*w)[0]*in[i-1] + (*w)[1]*in[i] + (*w)[2]*in[i+1];
174 } else { // mask it
175 out[i] = in[i];//use arbitrary value
176 }
177 outmask[i] = mask[i];
178 }
179}
180
181
182void mathutil::runningMedian(Vector<Float>& out, Vector<Bool>& outflag,
183 const Vector<Float>& in, const Vector<Bool>& flag,
184 float width)
185{
186 uInt hwidth = Int(width+0.5);
187 uInt fwidth = hwidth*2+1;
188 out.resize(in.nelements());
189 outflag.resize(flag.nelements());
190 MedianSlider ms(hwidth);
191 Slice sl(0, fwidth-1);
192 Float medval = ms.add(const_cast<Vector<Float>& >(in)(sl),
193 const_cast<Vector<Bool>& >(flag)(sl));
194 (void) medval;//suppress unused warning
195 uInt n = in.nelements();
196 for (uInt i=hwidth; i<(n-hwidth); ++i) {
197 // add data value
198 out[i] = ms.add(in[i+hwidth], flag[i+hwidth]);
199 outflag[i] = (ms.nval() > 0 );
200 }
201 // replicate edge values from first value with full width of values
202 for (uInt i=0;i<hwidth;++i) {
203 out[i] = out[hwidth];
204 outflag[i] = outflag[hwidth];
205 out[n-1-i] = out[n-1-hwidth];
206 outflag[n-1-i] = outflag[n-1-hwidth];
207 }
208}
209
210void mathutil::polyfit(Vector<Float>& out, Vector<Bool>& outmask,
211 const Vector<Float>& in, const Vector<Bool>& mask,
212 float width, int order)
213{
214 uInt hwidth = Int(width+0.5);
215 uInt fwidth = hwidth*2+1;
216 out.resize(in.nelements());
217 outmask.resize(mask.nelements());
218 LinearFit<Float> fitter;
219 Polynomial<Float> poly(order);
220 fitter.setFunction(poly);
221 Vector<Float> sigma(fwidth);
222 sigma = 1.0;
223 Vector<Float> parms;
224 Vector<Float> x(fwidth);
225 indgen(x);
226
227 uInt n = in.nelements();
228
229 for (uInt i=hwidth; i<(n-hwidth); ++i) {
230 // add data value
231 if (mask[i]) {
232 Slice sl(i-hwidth, fwidth);
233 const Vector<Float> &y = const_cast<Vector<Float>& >(in)(sl);
234 const Vector<Bool> &m = const_cast<Vector<Bool>& >(mask)(sl);
235 parms = fitter.fit(x, y, sigma, &m);
236
237 poly.setCoefficients(parms);
238 out[i] = poly(x[hwidth]);//cout << in[i] <<"->"<<out[i]<<endl;
239 } else {
240 out[i] = in[i];
241 }
242 outmask[i] = mask[i];
243 }
244 // replicate edge values from first value with full width of values
245 for (uInt i=0;i<hwidth;++i) {
246 out[i] = out[hwidth];
247 outmask[i] = outmask[hwidth];
248 out[n-1-i] = out[n-1-hwidth];
249 outmask[n-1-i] = outmask[n-1-hwidth];
250 }
251}
252
253void mathutil::doZeroOrderInterpolation(casa::Vector<casa::Float>& data,
254 std::vector<bool>& mask) {
255 int fstart = -1;
256 int fend = -1;
257 for (uInt i = 0; i < mask.size(); ++i) {
258 if (!mask[i]) {
259 fstart = i;
260 while (!mask[i] && i < mask.size()) {
261 fend = i;
262 i++;
263 }
264 }
265
266 // execute interpolation as the following criteria:
267 // (1) for a masked region inside the spectrum, replace the spectral
268 // values with the mean of those at the two channels just outside
269 // the both edges of the masked region.
270 // (2) for a masked region at the spectral edge, replace the values
271 // with the one at the nearest non-masked channel.
272 // (ZOH, but bilateral)
273 Float interp = 0.0;
274 if (fstart-1 > 0) {
275 interp = data[fstart-1];
276 if (fend+1 < Int(data.nelements())) {
277 interp = (interp + data[fend+1]) / 2.0;
278 }
279 } else {
280 interp = data[fend+1];
281 }
282 if (fstart > -1 && fend > -1) {
283 for (int j = fstart; j <= fend; ++j) {
284 data[j] = interp;
285 }
286 }
287
288 fstart = -1;
289 fend = -1;
290 }
291
292}
293
294void mathutil::rotateRA( Vector<Double> &v )
295{
296 uInt len = v.nelements() ;
297 Vector<Double> work( len ) ;
298
299 for ( uInt i = 0 ; i < len ; i++ ) {
300 work[i] = fmod( v[i], C::_2pi ) ;
301 if ( work[i] < 0.0 ) {
302 work[i] += C::_2pi ;
303 }
304 }
305 //cout << "zero2twopi: out=" << work << endl ;
306
307 Vector<uInt> quad( len ) ;
308 Vector<uInt> nquad( 4, 0 ) ;
309 for ( uInt i = 0 ; i < len ; i++ ) {
310 uInt q = uInt( work[i] / C::pi_2 ) ;
311 nquad[q]++ ;
312 quad[i] = q ;
313 }
314 //cout << "nquad = " << nquad << endl ;
315
316 Vector<Bool> rot( 4, False ) ;
317 if ( nquad[0] > 0 && nquad[3] > 0
318 && ( nquad[1] == 0 || nquad[2] == 0 ) ) {
319 //cout << "need rotation" << endl ;
320 rot[3] = True ;
321 rot[2] = (nquad[1]==0 && nquad[2]>0) ;
322 }
323 //cout << "rot=" << rot << endl ;
324
325 for ( uInt i = 0 ; i < len ; i++ ) {
326 if ( rot[quad[i]] ) {
327 v[i] = work[i] - C::_2pi ;
328 }
329 else {
330 v[i] = work[i] ;
331 }
332 }
333}
334
335void mathutil::rotateRA( const Vector<Double> &in,
336 Vector<Double> &out )
337{
338 out = in.copy() ;
339 rotateRA( out ) ;
340}
341
342double mathutil::gettimeofday_sec()
343{
344 struct timeval tv ;
345 gettimeofday( &tv, NULL ) ;
346 return tv.tv_sec + (double)tv.tv_usec*1.0e-6 ;
347}
Note: See TracBrowser for help on using the repository browser.