Line data Source code
1 : //# RFRowClipper.cc: this defines RFRowClipper
2 : //# Copyright (C) 2000,2001
3 : //# Associated Universities, Inc. Washington DC, USA.
4 : //#
5 : //# This library is free software; you can redistribute it and/or modify it
6 : //# under the terms of the GNU Library General Public License as published by
7 : //# the Free Software Foundation; either version 2 of the License, or (at your
8 : //# option) any later version.
9 : //#
10 : //# This library is distributed in the hope that it will be useful, but WITHOUT
11 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
13 : //# License for more details.
14 : //#
15 : //# You should have received a copy of the GNU Library General Public License
16 : //# along with this library; if not, write to the Free Software Foundation,
17 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
18 : //#
19 : //# Correspondence concerning AIPS++ should be addressed as follows:
20 : //# Internet email: aips2-request@nrao.edu.
21 : //# Postal address: AIPS++ Project Office
22 : //# National Radio Astronomy Observatory
23 : //# 520 Edgemont Road
24 : //# Charlottesville, VA 22903-2475 USA
25 : //#
26 : //# $Id$
27 : #include <flagging/Flagging/Flagger.h>
28 : #include <flagging/Flagging/RFFlagCube.h>
29 : #include <flagging/Flagging/RFChunkStats.h>
30 : #include <flagging/Flagging/RFRowClipper.h>
31 : #include <casacore/casa/Arrays/LogiVector.h>
32 : #include <casacore/casa/Arrays/ArrayMath.h>
33 : #include <casacore/casa/Arrays/MaskArrMath.h>
34 : #include <casacore/casa/Arrays/ArrayLogical.h>
35 : #include <casacore/casa/Arrays/Slice.h>
36 : #include <casacore/scimath/Mathematics/MedianSlider.h>
37 :
38 : using namespace casacore;
39 : namespace casa { //# NAMESPACE CASA - BEGIN
40 :
41 0 : RFRowClipper::RFRowClipper( RFChunkStats &ch,RFFlagCube &fl,Float clip,uInt hw,uInt maxp ) :
42 : chunk(ch),flag(fl),clip_level(clip),halfwin(hw),maxpass(maxp),
43 0 : os(fl.logSink())
44 0 : {}
45 :
46 0 : void RFRowClipper::init( uInt ni,uInt nt )
47 : {
48 0 : sig = Matrix<Float>(ntime=nt,nifr=ni,-1);
49 0 : sig0 = Matrix<Float>(nt,ni,-1);
50 0 : sigupdated = Vector<Bool>(ni,false);
51 0 : }
52 :
53 0 : void RFRowClipper::cleanup ()
54 : {
55 0 : sig.resize();
56 0 : sig0.resize();
57 0 : sigupdated.resize();
58 0 : }
59 :
60 0 : void RFRowClipper::reset ()
61 : {
62 0 : sigupdated = false;
63 0 : }
64 :
65 0 : Float RFRowClipper::updateSigma (uInt &ifrmax,uInt &itmax,Bool flag_rows, bool clear_flags )
66 : {
67 0 : Vector<Float> medsigma(ntime);
68 0 : Vector<Float> diffsigma(ntime);
69 0 : Vector<Float> diffs(ntime);
70 :
71 0 : Float dmax=0;
72 0 : ifrmax=itmax=0;
73 0 : RFlagWord fm = flag.flagMask()|flag.fullCorrMask();
74 :
75 0 : for( uInt ifr=0; ifr<nifr; ifr++ )
76 : {
77 0 : if( sigupdated(ifr) )
78 : {
79 : Bool fl;
80 : Float d;
81 0 : Vector<Float> sigma( sig.column(ifr) );
82 0 : Bool recalc=true;
83 0 : for( uInt ipass=0; ipass<maxpass && recalc; ipass++ ) // loop while some rows are being flagged
84 : {
85 0 : uInt idiff=0;
86 0 : recalc=false;
87 : // Precompute mask of valid sigmas: existing and not flagged
88 0 : LogicalVector valid(ntime,true);
89 0 : for( uInt i=0; i<ntime; i++ )
90 0 : if( sigma(i)<=0 || flag.getRowFlag(ifr,i)&fm )
91 0 : valid(i) = false;
92 :
93 : // If we have a valid half-window specified, then compute diff WRT
94 : // to a sliding median.
95 0 : if( halfwin>0 )
96 : {
97 0 : MedianSlider msl(halfwin);
98 0 : for( uInt it = 0; it<ntime; it++ )
99 : {
100 0 : msl.add( sigma(it), !valid(it) );
101 0 : if( it>=halfwin )
102 : {
103 0 : medsigma(it-halfwin) = msl.median();
104 0 : diffsigma(it-halfwin) = d = abs( msl.diff(fl) );
105 0 : if( !fl )
106 0 : diffs(idiff++) = d;
107 : }
108 : }
109 0 : for( uInt it=ntime-halfwin; it<ntime; it++ )
110 : {
111 0 : msl.next();
112 0 : medsigma(it) = msl.median();
113 0 : diffsigma(it) = d = abs(msl.diff(fl));
114 0 : if( !fl )
115 0 : diffs(idiff++) = d;
116 : }
117 : }
118 : else // No half-window, compute diff WRT global median
119 : {
120 0 : Vector<Float> s;
121 0 : s = sigma(valid);
122 0 : Float med = median(s);
123 0 : medsigma.set(med);
124 0 : diffsigma = abs( sigma - med );
125 0 : diffs.resize();
126 0 : diffs = diffsigma(valid);
127 0 : idiff = diffs.nelements();
128 : }
129 0 : if( !idiff ) // no data? go on
130 0 : continue;
131 : // compute threshold, using median of the good datums
132 0 : Float meddiff = idiff ? median( diffs( Slice(0,idiff) ) ) : 0;
133 0 : Float thr = clip_level*meddiff;
134 0 : uInt nbad=0;
135 0 : LogicalVector goodsigma( diffsigma<thr );
136 0 : for( uInt it=0; it<ntime; it++ )
137 : {
138 0 : Float s=sigma(it);
139 0 : if( s>0 )
140 : {
141 : // for good rows (or when not using row flagging at all)
142 : // update stats and clear flags, if needed
143 0 : if( !flag_rows || goodsigma(it) )
144 : {
145 0 : Bool res = false;
146 0 : if( flag_rows && clear_flags ) // clear row flag
147 : {
148 0 : recalc |= ( res = flag.clearRowFlag(ifr,it) );
149 0 : for( uInt ich=0; ich<chunk.num(CHAN); ich++ )
150 0 : flag.clearFlag(ich,ifr);
151 : }
152 0 : Float s0 = sig0(it,ifr),
153 0 : m = max(s,s0),
154 0 : d= m!=0 ? abs(s-s0)/m : 0;
155 0 : if( d>dmax ) // compute new max difference in sigma
156 : {
157 0 : dmax=d;
158 0 : ifrmax=ifr; itmax=it;
159 : }
160 : }
161 : else // set flags on apparently bad rows
162 : {
163 0 : Bool res = flag.setRowFlag(ifr,it);
164 0 : recalc |= res;
165 0 : for( uInt ich=0; ich<chunk.num(CHAN); ich++ )
166 0 : flag.setFlag(ich,ifr);
167 0 : nbad++;
168 : }
169 : }
170 : else // ignore rows that are apriori bad/nonexistent
171 : {
172 : }
173 : }
174 0 : String ifrid( chunk.ifrString(ifr) );
175 : // dprintf(os,"IFR %d (%s): %d rows flagged, recalc=%d\n",ifr,ifrid.chars(),nbad,(Int)recalc);
176 : } // endwhile(recalc)
177 : } // endif( sigupated(ifr) )
178 : } // endfor( ifr )
179 :
180 0 : sig0 = sig;
181 :
182 : // dprintf(os,"Max diff (%f) at ifr %d (%s), it %d: new sigma is %f\n",
183 : // dmax,ifrmax,chunk.ifrString(ifrmax).chars(),itmax,sig0(itmax,ifrmax));
184 :
185 0 : return dmax;
186 : }
187 :
188 : } //# NAMESPACE CASA - END
189 :
|
