[2] | 1 | //#---------------------------------------------------------------------------
|
---|
| 2 | //# SDMath.cc: A collection of single dish mathematical operations
|
---|
| 3 | //#---------------------------------------------------------------------------
|
---|
| 4 | //# Copyright (C) 2004
|
---|
[125] | 5 | //# ATNF
|
---|
[2] | 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 | //#---------------------------------------------------------------------------
|
---|
[38] | 31 | #include <vector>
|
---|
| 32 |
|
---|
[81] | 33 | #include <casa/aips.h>
|
---|
| 34 | #include <casa/BasicSL/String.h>
|
---|
| 35 | #include <casa/Arrays/IPosition.h>
|
---|
| 36 | #include <casa/Arrays/Array.h>
|
---|
[130] | 37 | #include <casa/Arrays/ArrayIter.h>
|
---|
| 38 | #include <casa/Arrays/VectorIter.h>
|
---|
[81] | 39 | #include <casa/Arrays/ArrayMath.h>
|
---|
| 40 | #include <casa/Arrays/ArrayLogical.h>
|
---|
| 41 | #include <casa/Arrays/MaskedArray.h>
|
---|
| 42 | #include <casa/Arrays/MaskArrMath.h>
|
---|
| 43 | #include <casa/Arrays/MaskArrLogi.h>
|
---|
[130] | 44 | #include <casa/Exceptions.h>
|
---|
[2] | 45 |
|
---|
[81] | 46 | #include <tables/Tables/Table.h>
|
---|
| 47 | #include <tables/Tables/ScalarColumn.h>
|
---|
| 48 | #include <tables/Tables/ArrayColumn.h>
|
---|
[2] | 49 |
|
---|
[130] | 50 | #include <lattices/Lattices/LatticeUtilities.h>
|
---|
| 51 | #include <lattices/Lattices/RebinLattice.h>
|
---|
[81] | 52 | #include <coordinates/Coordinates/SpectralCoordinate.h>
|
---|
[130] | 53 | #include <coordinates/Coordinates/CoordinateSystem.h>
|
---|
| 54 | #include <coordinates/Coordinates/CoordinateUtil.h>
|
---|
[38] | 55 |
|
---|
| 56 | #include "MathUtils.h"
|
---|
[2] | 57 | #include "SDContainer.h"
|
---|
| 58 | #include "SDMemTable.h"
|
---|
| 59 |
|
---|
| 60 | #include "SDMath.h"
|
---|
| 61 |
|
---|
[125] | 62 | using namespace casa;
|
---|
[83] | 63 | using namespace asap;
|
---|
[2] | 64 |
|
---|
[144] | 65 | CountedPtr<SDMemTable> SDMath::average (const Block<CountedPtr<SDMemTable> >& in,
|
---|
| 66 | const Vector<Bool>& mask, bool scanAv,
|
---|
| 67 | const std::string& weightStr)
|
---|
[130] | 68 | //
|
---|
[144] | 69 | // Weighted averaging of spectra from one or more Tables.
|
---|
[130] | 70 | //
|
---|
| 71 | {
|
---|
[2] | 72 |
|
---|
[163] | 73 | // Convert weight type
|
---|
| 74 |
|
---|
| 75 | WeightType wtType = NONE;
|
---|
| 76 | convertWeightString (wtType, weightStr);
|
---|
| 77 |
|
---|
[144] | 78 | // Create output Table by cloning from the first table
|
---|
[2] | 79 |
|
---|
[144] | 80 | SDMemTable* pTabOut = new SDMemTable(*in[0],True);
|
---|
[130] | 81 |
|
---|
[144] | 82 | // Setup
|
---|
[130] | 83 |
|
---|
[144] | 84 | const uInt axis = 3; // Spectral axis
|
---|
| 85 | IPosition shp = in[0]->rowAsMaskedArray(0).shape(); // Must not change
|
---|
| 86 | Array<Float> arr(shp);
|
---|
| 87 | Array<Bool> barr(shp);
|
---|
| 88 | const Bool useMask = (mask.nelements() == shp(axis));
|
---|
[130] | 89 |
|
---|
[144] | 90 | // Columns from Tables
|
---|
[130] | 91 |
|
---|
[144] | 92 | ROArrayColumn<Float> tSysCol;
|
---|
| 93 | ROScalarColumn<Double> mjdCol;
|
---|
| 94 | ROScalarColumn<String> srcNameCol;
|
---|
| 95 | ROScalarColumn<Double> intCol;
|
---|
| 96 | ROArrayColumn<uInt> fqIDCol;
|
---|
[130] | 97 |
|
---|
[144] | 98 | // Create accumulation MaskedArray. We accumulate for each channel,if,pol,beam
|
---|
| 99 | // Note that the mask of the accumulation array will ALWAYS remain ALL True.
|
---|
| 100 | // The MA is only used so that when data which is masked Bad is added to it,
|
---|
| 101 | // that data does not contribute.
|
---|
| 102 |
|
---|
| 103 | Array<Float> zero(shp);
|
---|
| 104 | zero=0.0;
|
---|
| 105 | Array<Bool> good(shp);
|
---|
| 106 | good = True;
|
---|
| 107 | MaskedArray<Float> sum(zero,good);
|
---|
| 108 |
|
---|
| 109 | // Counter arrays
|
---|
| 110 |
|
---|
| 111 | Array<Float> nPts(shp); // Number of points
|
---|
| 112 | nPts = 0.0;
|
---|
| 113 | Array<Float> nInc(shp); // Increment
|
---|
| 114 | nInc = 1.0;
|
---|
| 115 |
|
---|
| 116 | // Create accumulation Array for variance. We accumulate for
|
---|
| 117 | // each if,pol,beam, but average over channel. So we need
|
---|
| 118 | // a shape with one less axis dropping channels.
|
---|
| 119 |
|
---|
| 120 | const uInt nAxesSub = shp.nelements() - 1;
|
---|
| 121 | IPosition shp2(nAxesSub);
|
---|
| 122 | for (uInt i=0,j=0; i<(nAxesSub+1); i++) {
|
---|
| 123 | if (i!=axis) {
|
---|
| 124 | shp2(j) = shp(i);
|
---|
| 125 | j++;
|
---|
| 126 | }
|
---|
[2] | 127 | }
|
---|
[144] | 128 | Array<Float> sumSq(shp2);
|
---|
| 129 | sumSq = 0.0;
|
---|
| 130 | IPosition pos2(nAxesSub,0); // For indexing
|
---|
[130] | 131 |
|
---|
[144] | 132 | // Time-related accumulators
|
---|
[130] | 133 |
|
---|
[144] | 134 | Double time;
|
---|
| 135 | Double timeSum = 0.0;
|
---|
| 136 | Double intSum = 0.0;
|
---|
| 137 | Double interval = 0.0;
|
---|
[130] | 138 |
|
---|
[144] | 139 | // To get the right shape for the Tsys accumulator we need to
|
---|
| 140 | // access a column from the first table. The shape of this
|
---|
| 141 | // array must not change
|
---|
[130] | 142 |
|
---|
[144] | 143 | Array<Float> tSysSum;
|
---|
| 144 | {
|
---|
| 145 | const Table& tabIn = in[0]->table();
|
---|
| 146 | tSysCol.attach(tabIn,"TSYS");
|
---|
| 147 | tSysSum.resize(tSysCol.shape(0));
|
---|
| 148 | }
|
---|
| 149 | tSysSum =0.0;
|
---|
| 150 | Array<Float> tSys;
|
---|
| 151 |
|
---|
| 152 | // Scan and row tracking
|
---|
| 153 |
|
---|
| 154 | Int oldScanID = 0;
|
---|
| 155 | Int outScanID = 0;
|
---|
| 156 | Int scanID = 0;
|
---|
| 157 | Int rowStart = 0;
|
---|
| 158 | Int nAccum = 0;
|
---|
| 159 | Int tableStart = 0;
|
---|
| 160 |
|
---|
| 161 | // Source and FreqID
|
---|
| 162 |
|
---|
| 163 | String sourceName, oldSourceName, sourceNameStart;
|
---|
| 164 | Vector<uInt> freqID, freqIDStart, oldFreqID;
|
---|
| 165 |
|
---|
| 166 | // Loop over tables
|
---|
| 167 |
|
---|
| 168 | Float fac = 1.0;
|
---|
| 169 | const uInt nTables = in.nelements();
|
---|
| 170 | for (uInt iTab=0; iTab<nTables; iTab++) {
|
---|
| 171 |
|
---|
| 172 | // Attach columns to Table
|
---|
| 173 |
|
---|
| 174 | const Table& tabIn = in[iTab]->table();
|
---|
| 175 | tSysCol.attach(tabIn, "TSYS");
|
---|
| 176 | mjdCol.attach(tabIn, "TIME");
|
---|
| 177 | srcNameCol.attach(tabIn, "SRCNAME");
|
---|
| 178 | intCol.attach(tabIn, "INTERVAL");
|
---|
| 179 | fqIDCol.attach(tabIn, "FREQID");
|
---|
| 180 |
|
---|
| 181 | // Loop over rows in Table
|
---|
| 182 |
|
---|
| 183 | const uInt nRows = in[iTab]->nRow();
|
---|
| 184 | for (uInt iRow=0; iRow<nRows; iRow++) {
|
---|
| 185 |
|
---|
| 186 | // Check conformance
|
---|
| 187 |
|
---|
| 188 | IPosition shp2 = in[iTab]->rowAsMaskedArray(iRow).shape();
|
---|
| 189 | if (!shp.isEqual(shp2)) {
|
---|
| 190 | throw (AipsError("Shapes for all rows must be the same"));
|
---|
| 191 | }
|
---|
| 192 |
|
---|
| 193 | // If we are not doing scan averages, make checks for source and
|
---|
| 194 | // frequency setup and warn if averaging across them
|
---|
| 195 |
|
---|
| 196 | // Get copy of Scan Container for this row
|
---|
| 197 |
|
---|
| 198 | SDContainer sc = in[iTab]->getSDContainer(iRow);
|
---|
| 199 | scanID = sc.scanid;
|
---|
| 200 |
|
---|
| 201 | // Get quantities from columns
|
---|
| 202 |
|
---|
| 203 | srcNameCol.getScalar(iRow, sourceName);
|
---|
| 204 | mjdCol.get(iRow, time);
|
---|
| 205 | tSysCol.get(iRow, tSys);
|
---|
| 206 | intCol.get(iRow, interval);
|
---|
| 207 | fqIDCol.get(iRow, freqID);
|
---|
| 208 |
|
---|
| 209 | // Initialize first source and freqID
|
---|
| 210 |
|
---|
| 211 | if (iRow==0 && iTab==0) {
|
---|
| 212 | sourceNameStart = sourceName;
|
---|
| 213 | freqIDStart = freqID;
|
---|
| 214 | }
|
---|
| 215 |
|
---|
| 216 | // If we are doing scan averages, see if we are at the end of an
|
---|
| 217 | // accumulation period (scan). We must check soutce names too,
|
---|
| 218 | // since we might have two tables with one scan each but different
|
---|
| 219 | // source names; we shouldn't average different sources together
|
---|
| 220 |
|
---|
| 221 | if (scanAv && ( (scanID != oldScanID) ||
|
---|
| 222 | (iRow==0 && iTab>0 && sourceName!=oldSourceName))) {
|
---|
| 223 |
|
---|
| 224 | // Normalize data in 'sum' accumulation array according to weighting scheme
|
---|
| 225 |
|
---|
| 226 | normalize (sum, sumSq, nPts, wtType, axis, nAxesSub);
|
---|
| 227 |
|
---|
| 228 | // Fill scan container. The source and freqID come from the
|
---|
| 229 | // first row of the first table that went into this average (
|
---|
| 230 | // should be the same for all rows in the scan average)
|
---|
| 231 |
|
---|
| 232 | Float nR(nAccum);
|
---|
| 233 | fillSDC (sc, sum.getMask(), sum.getArray(), tSysSum/nR, outScanID,
|
---|
| 234 | timeSum/nR, intSum, sourceNameStart, freqIDStart);
|
---|
| 235 |
|
---|
| 236 | // Write container out to Table
|
---|
| 237 |
|
---|
| 238 | pTabOut->putSDContainer(sc);
|
---|
| 239 |
|
---|
| 240 | // Reset accumulators
|
---|
| 241 |
|
---|
| 242 | sum = 0.0;
|
---|
| 243 | sumSq = 0.0;
|
---|
| 244 | nAccum = 0;
|
---|
| 245 | //
|
---|
| 246 | tSysSum =0.0;
|
---|
| 247 | timeSum = 0.0;
|
---|
| 248 | intSum = 0.0;
|
---|
| 249 |
|
---|
| 250 | // Increment
|
---|
| 251 |
|
---|
| 252 | rowStart = iRow; // First row for next accumulation
|
---|
| 253 | tableStart = iTab; // First table for next accumulation
|
---|
| 254 | sourceNameStart = sourceName; // First source name for next accumulation
|
---|
| 255 | freqIDStart = freqID; // First FreqID for next accumulation
|
---|
| 256 | //
|
---|
| 257 | oldScanID = scanID;
|
---|
| 258 | outScanID += 1; // Scan ID for next accumulation period
|
---|
| 259 | }
|
---|
| 260 |
|
---|
[146] | 261 | // Accumulate
|
---|
[144] | 262 |
|
---|
[146] | 263 | accumulate (timeSum, intSum, nAccum, sum, sumSq, nPts, tSysSum,
|
---|
| 264 | tSys, nInc, mask, time, interval, in, iTab, iRow, axis,
|
---|
| 265 | nAxesSub, useMask, wtType);
|
---|
[144] | 266 | //
|
---|
| 267 | oldSourceName = sourceName;
|
---|
| 268 | oldFreqID = freqID;
|
---|
| 269 | }
|
---|
| 270 | }
|
---|
| 271 |
|
---|
| 272 | // OK at this point we have accumulation data which is either
|
---|
| 273 | // - accumulated from all tables into one row
|
---|
| 274 | // or
|
---|
| 275 | // - accumulated from the last scan average
|
---|
| 276 | //
|
---|
| 277 | // Normalize data in 'sum' accumulation array according to weighting scheme
|
---|
| 278 |
|
---|
| 279 | normalize (sum, sumSq, nPts, wtType, axis, nAxesSub);
|
---|
| 280 |
|
---|
| 281 | // Create and fill container. The container we clone will be from
|
---|
| 282 | // the last Table and the first row that went into the current
|
---|
| 283 | // accumulation. It probably doesn't matter that much really...
|
---|
| 284 |
|
---|
| 285 | Float nR(nAccum);
|
---|
| 286 | SDContainer sc = in[tableStart]->getSDContainer(rowStart);
|
---|
| 287 | fillSDC (sc, sum.getMask(), sum.getArray(), tSysSum/nR, outScanID,
|
---|
| 288 | timeSum/nR, intSum, sourceNameStart, freqIDStart);
|
---|
| 289 | //
|
---|
| 290 | pTabOut->putSDContainer(sc);
|
---|
| 291 | /*
|
---|
| 292 | cout << endl;
|
---|
| 293 | cout << "Last accumulation for output scan ID " << outScanID << endl;
|
---|
| 294 | cout << " The first row in this accumulation is " << rowStart << endl;
|
---|
| 295 | cout << " The number of rows accumulated is " << nAccum << endl;
|
---|
| 296 | cout << " The first table in this accumulation is " << tableStart << endl;
|
---|
| 297 | cout << " The first source in this accumulation is " << sourceNameStart << endl;
|
---|
| 298 | cout << " The first freqID in this accumulation is " << freqIDStart << endl;
|
---|
| 299 | cout << " Average time stamp = " << timeSum/nR << endl;
|
---|
| 300 | cout << " Integrated time = " << intSum << endl;
|
---|
| 301 | */
|
---|
| 302 | return CountedPtr<SDMemTable>(pTabOut);
|
---|
[2] | 303 | }
|
---|
[9] | 304 |
|
---|
[144] | 305 |
|
---|
| 306 |
|
---|
[85] | 307 | CountedPtr<SDMemTable>
|
---|
| 308 | SDMath::quotient(const CountedPtr<SDMemTable>& on,
|
---|
[130] | 309 | const CountedPtr<SDMemTable>& off)
|
---|
| 310 | //
|
---|
| 311 | // Compute quotient spectrum
|
---|
| 312 | //
|
---|
| 313 | {
|
---|
| 314 | const uInt nRows = on->nRow();
|
---|
| 315 | if (off->nRow() != nRows) {
|
---|
| 316 | throw (AipsError("Input Scan Tables must have the same number of rows"));
|
---|
| 317 | }
|
---|
[85] | 318 |
|
---|
[130] | 319 | // Input Tables and columns
|
---|
| 320 |
|
---|
[9] | 321 | Table ton = on->table();
|
---|
| 322 | Table toff = off->table();
|
---|
[85] | 323 | ROArrayColumn<Float> tsys(toff, "TSYS");
|
---|
[9] | 324 | ROScalarColumn<Double> mjd(ton, "TIME");
|
---|
[15] | 325 | ROScalarColumn<Double> integr(ton, "INTERVAL");
|
---|
[9] | 326 | ROScalarColumn<String> srcn(ton, "SRCNAME");
|
---|
[38] | 327 | ROArrayColumn<uInt> freqidc(ton, "FREQID");
|
---|
| 328 |
|
---|
[130] | 329 | // Output Table cloned from input
|
---|
[85] | 330 |
|
---|
[15] | 331 | SDMemTable* sdmt = new SDMemTable(*on, True);
|
---|
[130] | 332 |
|
---|
| 333 | // Loop over rows
|
---|
| 334 |
|
---|
| 335 | for (uInt i=0; i<nRows; i++) {
|
---|
| 336 | MaskedArray<Float> mon(on->rowAsMaskedArray(i));
|
---|
| 337 | MaskedArray<Float> moff(off->rowAsMaskedArray(i));
|
---|
| 338 | IPosition ipon = mon.shape();
|
---|
| 339 | IPosition ipoff = moff.shape();
|
---|
| 340 | //
|
---|
| 341 | Array<Float> tsarr;
|
---|
| 342 | tsys.get(i, tsarr);
|
---|
| 343 | if (ipon != ipoff && ipon != tsarr.shape()) {
|
---|
| 344 | throw(AipsError("on/off not conformant"));
|
---|
| 345 | }
|
---|
| 346 |
|
---|
| 347 | // Compute quotient
|
---|
| 348 |
|
---|
| 349 | MaskedArray<Float> tmp = (mon-moff);
|
---|
| 350 | Array<Float> out(tmp.getArray());
|
---|
| 351 | out /= moff;
|
---|
| 352 | out *= tsarr;
|
---|
[163] | 353 | Array<Bool> outflagsb = mon.getMask() && moff.getMask();
|
---|
[130] | 354 |
|
---|
| 355 | // Fill container for this row
|
---|
| 356 |
|
---|
| 357 | SDContainer sc = on->getSDContainer();
|
---|
[163] | 358 | //
|
---|
| 359 | putDataInSDC (sc, out, outflagsb);
|
---|
[130] | 360 | sc.putTsys(tsarr);
|
---|
| 361 | sc.scanid = 0;
|
---|
| 362 |
|
---|
| 363 | // Put new row in output Table
|
---|
| 364 |
|
---|
| 365 | sdmt->putSDContainer(sc);
|
---|
| 366 | }
|
---|
| 367 | //
|
---|
[9] | 368 | return CountedPtr<SDMemTable>(sdmt);
|
---|
| 369 | }
|
---|
[48] | 370 |
|
---|
[152] | 371 | void SDMath::multiplyInSitu(SDMemTable* pIn, Float factor, Bool doAll)
|
---|
[146] | 372 | {
|
---|
[152] | 373 | const uInt what = 0;
|
---|
| 374 | SDMemTable* pOut = localOperate (*pIn, factor, doAll, what);
|
---|
[146] | 375 | *pIn = *pOut;
|
---|
| 376 | delete pOut;
|
---|
| 377 | }
|
---|
| 378 |
|
---|
| 379 |
|
---|
[85] | 380 | CountedPtr<SDMemTable>
|
---|
[152] | 381 | SDMath::multiply(const CountedPtr<SDMemTable>& in, Float factor, Bool doAll)
|
---|
[130] | 382 | {
|
---|
[152] | 383 | const uInt what = 0;
|
---|
| 384 | return CountedPtr<SDMemTable>(localOperate (*in, factor, doAll, what));
|
---|
[15] | 385 | }
|
---|
[48] | 386 |
|
---|
[152] | 387 |
|
---|
| 388 | void SDMath::addInSitu (SDMemTable* pIn, Float offset, Bool doAll)
|
---|
| 389 | {
|
---|
| 390 | const uInt what = 1;
|
---|
| 391 | SDMemTable* pOut = localOperate (*pIn, offset, doAll, what);
|
---|
| 392 | *pIn = *pOut;
|
---|
| 393 | delete pOut;
|
---|
| 394 | }
|
---|
| 395 |
|
---|
| 396 |
|
---|
[107] | 397 | CountedPtr<SDMemTable>
|
---|
[152] | 398 | SDMath::add(const CountedPtr<SDMemTable>& in, Float offset, Bool doAll)
|
---|
[130] | 399 | {
|
---|
[152] | 400 | const uInt what = 1;
|
---|
| 401 | return CountedPtr<SDMemTable>(localOperate(*in, offset, doAll, what));
|
---|
[107] | 402 | }
|
---|
| 403 |
|
---|
| 404 |
|
---|
[130] | 405 | CountedPtr<SDMemTable>
|
---|
| 406 | SDMath::hanning(const CountedPtr<SDMemTable>& in)
|
---|
| 407 | //
|
---|
| 408 | // Hanning smooth each row
|
---|
| 409 | // Should Tsys be smoothed ?
|
---|
| 410 | //
|
---|
| 411 | {
|
---|
| 412 | SDMemTable* sdmt = new SDMemTable(*in,True);
|
---|
[107] | 413 |
|
---|
[130] | 414 | // Loop over rows in Table
|
---|
[48] | 415 |
|
---|
[130] | 416 | for (uInt ri=0; ri < in->nRow(); ++ri) {
|
---|
[38] | 417 |
|
---|
[130] | 418 | // Get data
|
---|
[125] | 419 |
|
---|
[130] | 420 | const MaskedArray<Float>& marr(in->rowAsMaskedArray(ri));
|
---|
[85] | 421 | Array<Float> arr = marr.getArray();
|
---|
| 422 | Array<Bool> barr = marr.getMask();
|
---|
[130] | 423 |
|
---|
| 424 | // Smooth along the channels axis
|
---|
| 425 |
|
---|
| 426 | uInt axis = 3;
|
---|
| 427 | VectorIterator<Float> itData(arr, axis);
|
---|
| 428 | VectorIterator<Bool> itMask(barr, axis);
|
---|
| 429 | Vector<Float> outv;
|
---|
| 430 | Vector<Bool> outm;
|
---|
| 431 | while (!itData.pastEnd()) {
|
---|
| 432 | mathutil::hanning(outv, outm, itData.vector(), itMask.vector());
|
---|
| 433 | itData.vector() = outv;
|
---|
| 434 | itMask.vector() = outm;
|
---|
| 435 | //
|
---|
| 436 | itData.next();
|
---|
| 437 | itMask.next();
|
---|
[38] | 438 | }
|
---|
[130] | 439 |
|
---|
| 440 | // Create and put back
|
---|
| 441 |
|
---|
[85] | 442 | SDContainer sc = in->getSDContainer(ri);
|
---|
[163] | 443 | putDataInSDC (sc, arr, barr);
|
---|
| 444 | //
|
---|
[85] | 445 | sdmt->putSDContainer(sc);
|
---|
[38] | 446 | }
|
---|
| 447 | return CountedPtr<SDMemTable>(sdmt);
|
---|
| 448 | }
|
---|
| 449 |
|
---|
[165] | 450 | void SDMath::averagePolInSitu (SDMemTable* pIn, const Vector<Bool>& mask)
|
---|
[130] | 451 | {
|
---|
[165] | 452 | SDMemTable* pOut = localAveragePol (*pIn, mask);
|
---|
| 453 | *pIn = *pOut;
|
---|
| 454 | delete pOut;
|
---|
| 455 | }
|
---|
[163] | 456 |
|
---|
[130] | 457 |
|
---|
[165] | 458 | CountedPtr<SDMemTable> SDMath::averagePol (const CountedPtr<SDMemTable>& in,
|
---|
| 459 | const Vector<Bool>& mask)
|
---|
| 460 | {
|
---|
| 461 | return CountedPtr<SDMemTable>(localAveragePol(*in, mask));
|
---|
[48] | 462 | }
|
---|
| 463 |
|
---|
[167] | 464 |
|
---|
| 465 | void SDMath::binInSitu (SDMemTable* pIn, int width)
|
---|
[130] | 466 | {
|
---|
[167] | 467 | const uInt what = 1;
|
---|
| 468 | SDMemTable* pOut = localBin (*pIn, Int(width));
|
---|
| 469 | *pIn = *pOut;
|
---|
| 470 | delete pOut;
|
---|
| 471 | }
|
---|
[85] | 472 |
|
---|
[130] | 473 |
|
---|
[167] | 474 | CountedPtr<SDMemTable> SDMath::bin (const CountedPtr<SDMemTable>& in, int width)
|
---|
| 475 | {
|
---|
| 476 | const uInt what = 1;
|
---|
| 477 | return CountedPtr<SDMemTable>(localBin(*in, Int(width)));
|
---|
| 478 | }
|
---|
[130] | 479 |
|
---|
| 480 |
|
---|
[85] | 481 |
|
---|
[130] | 482 |
|
---|
| 483 | std::vector<float> SDMath::statistic (const CountedPtr<SDMemTable>& in,
|
---|
| 484 | const std::vector<bool>& mask,
|
---|
| 485 | const std::string& which)
|
---|
| 486 | //
|
---|
| 487 | // Perhaps iteration over pol/beam/if should be in here
|
---|
| 488 | // and inside the nrow iteration ?
|
---|
| 489 | //
|
---|
| 490 | {
|
---|
| 491 | const uInt nRow = in->nRow();
|
---|
| 492 | std::vector<float> result(nRow);
|
---|
| 493 | Vector<Bool> msk(mask);
|
---|
| 494 |
|
---|
| 495 | // Specify cursor location
|
---|
| 496 |
|
---|
[152] | 497 | IPosition start, end;
|
---|
| 498 | getCursorLocation (start, end, *in);
|
---|
[130] | 499 |
|
---|
| 500 | // Loop over rows
|
---|
| 501 |
|
---|
| 502 | const uInt nEl = msk.nelements();
|
---|
| 503 | for (uInt ii=0; ii < in->nRow(); ++ii) {
|
---|
| 504 |
|
---|
| 505 | // Get row and deconstruct
|
---|
| 506 |
|
---|
| 507 | MaskedArray<Float> marr(in->rowAsMaskedArray(ii));
|
---|
| 508 | Array<Float> arr = marr.getArray();
|
---|
| 509 | Array<Bool> barr = marr.getMask();
|
---|
| 510 |
|
---|
| 511 | // Access desired piece of data
|
---|
| 512 |
|
---|
| 513 | Array<Float> v((arr(start,end)).nonDegenerate());
|
---|
| 514 | Array<Bool> m((barr(start,end)).nonDegenerate());
|
---|
| 515 |
|
---|
| 516 | // Apply OTF mask
|
---|
| 517 |
|
---|
| 518 | MaskedArray<Float> tmp;
|
---|
| 519 | if (m.nelements()==nEl) {
|
---|
| 520 | tmp.setData(v,m&&msk);
|
---|
| 521 | } else {
|
---|
| 522 | tmp.setData(v,m);
|
---|
| 523 | }
|
---|
| 524 |
|
---|
| 525 | // Get statistic
|
---|
| 526 |
|
---|
[144] | 527 | result[ii] = mathutil::statistics(which, tmp);
|
---|
[130] | 528 | }
|
---|
| 529 | //
|
---|
| 530 | return result;
|
---|
| 531 | }
|
---|
| 532 |
|
---|
[146] | 533 |
|
---|
| 534 |
|
---|
| 535 | // 'private' functions
|
---|
| 536 |
|
---|
[144] | 537 | void SDMath::fillSDC (SDContainer& sc,
|
---|
| 538 | const Array<Bool>& mask,
|
---|
| 539 | const Array<Float>& data,
|
---|
| 540 | const Array<Float>& tSys,
|
---|
| 541 | Int scanID, Double timeStamp,
|
---|
| 542 | Double interval, const String& sourceName,
|
---|
| 543 | const Vector<uInt>& freqID)
|
---|
| 544 | {
|
---|
[163] | 545 | // Data and mask
|
---|
| 546 |
|
---|
| 547 | putDataInSDC (sc, data, mask);
|
---|
| 548 |
|
---|
| 549 | // TSys
|
---|
| 550 |
|
---|
[144] | 551 | sc.putTsys(tSys);
|
---|
[130] | 552 |
|
---|
[144] | 553 | // Time things
|
---|
| 554 |
|
---|
| 555 | sc.timestamp = timeStamp;
|
---|
| 556 | sc.interval = interval;
|
---|
| 557 | sc.scanid = scanID;
|
---|
| 558 | //
|
---|
| 559 | sc.sourcename = sourceName;
|
---|
| 560 | sc.putFreqMap(freqID);
|
---|
| 561 | }
|
---|
| 562 |
|
---|
| 563 | void SDMath::normalize (MaskedArray<Float>& sum,
|
---|
| 564 | const Array<Float>& sumSq,
|
---|
| 565 | const Array<Float>& nPts,
|
---|
[163] | 566 | WeightType wtType, Int axis,
|
---|
[144] | 567 | Int nAxesSub)
|
---|
[130] | 568 | {
|
---|
[144] | 569 | IPosition pos2(nAxesSub,0);
|
---|
| 570 | //
|
---|
| 571 | if (wtType==NONE) {
|
---|
| 572 |
|
---|
| 573 | // We just average by the number of points accumulated.
|
---|
| 574 | // We need to make a MA out of nPts so that no divide by
|
---|
| 575 | // zeros occur
|
---|
| 576 |
|
---|
| 577 | MaskedArray<Float> t(nPts, (nPts>Float(0.0)));
|
---|
| 578 | sum /= t;
|
---|
| 579 | } else if (wtType==VAR) {
|
---|
| 580 |
|
---|
| 581 | // Normalize each spectrum by sum(1/var) where the variance
|
---|
| 582 | // is worked out for each spectrum
|
---|
| 583 |
|
---|
| 584 | Array<Float>& data = sum.getRWArray();
|
---|
| 585 | VectorIterator<Float> itData(data, axis);
|
---|
| 586 | while (!itData.pastEnd()) {
|
---|
| 587 | pos2 = itData.pos().getFirst(nAxesSub);
|
---|
| 588 | itData.vector() /= sumSq(pos2);
|
---|
| 589 | itData.next();
|
---|
| 590 | }
|
---|
| 591 | } else if (wtType==TSYS) {
|
---|
[130] | 592 | }
|
---|
| 593 | }
|
---|
[144] | 594 |
|
---|
[146] | 595 |
|
---|
| 596 | void SDMath::accumulate (Double& timeSum, Double& intSum, Int& nAccum,
|
---|
| 597 | MaskedArray<Float>& sum, Array<Float>& sumSq,
|
---|
| 598 | Array<Float>& nPts, Array<Float>& tSysSum,
|
---|
| 599 | const Array<Float>& tSys, const Array<Float>& nInc,
|
---|
| 600 | const Vector<Bool>& mask, Double time, Double interval,
|
---|
| 601 | const Block<CountedPtr<SDMemTable> >& in,
|
---|
| 602 | uInt iTab, uInt iRow, uInt axis,
|
---|
| 603 | uInt nAxesSub, Bool useMask,
|
---|
[163] | 604 | WeightType wtType)
|
---|
[146] | 605 | {
|
---|
| 606 |
|
---|
| 607 | // Get data
|
---|
| 608 |
|
---|
| 609 | MaskedArray<Float> dataIn(in[iTab]->rowAsMaskedArray(iRow));
|
---|
| 610 | Array<Float>& valuesIn = dataIn.getRWArray(); // writable reference
|
---|
| 611 | const Array<Bool>& maskIn = dataIn.getMask(); // RO reference
|
---|
| 612 | //
|
---|
| 613 | if (wtType==NONE) {
|
---|
| 614 | const MaskedArray<Float> n(nInc,dataIn.getMask());
|
---|
| 615 | nPts += n; // Only accumulates where mask==T
|
---|
| 616 | } else if (wtType==VAR) {
|
---|
| 617 |
|
---|
| 618 | // We are going to average the data, weighted by the noise for each pol, beam and IF.
|
---|
| 619 | // So therefore we need to iterate through by spectrum (axis 3)
|
---|
| 620 |
|
---|
| 621 | VectorIterator<Float> itData(valuesIn, axis);
|
---|
| 622 | ReadOnlyVectorIterator<Bool> itMask(maskIn, axis);
|
---|
| 623 | Float fac = 1.0;
|
---|
| 624 | IPosition pos(nAxesSub,0);
|
---|
| 625 | //
|
---|
| 626 | while (!itData.pastEnd()) {
|
---|
| 627 |
|
---|
| 628 | // Make MaskedArray of Vector, optionally apply OTF mask, and find scaling factor
|
---|
| 629 |
|
---|
| 630 | if (useMask) {
|
---|
| 631 | MaskedArray<Float> tmp(itData.vector(),mask&&itMask.vector());
|
---|
| 632 | fac = 1.0/variance(tmp);
|
---|
| 633 | } else {
|
---|
| 634 | MaskedArray<Float> tmp(itData.vector(),itMask.vector());
|
---|
| 635 | fac = 1.0/variance(tmp);
|
---|
| 636 | }
|
---|
| 637 |
|
---|
| 638 | // Scale data
|
---|
| 639 |
|
---|
| 640 | itData.vector() *= fac; // Writes back into 'dataIn'
|
---|
| 641 | //
|
---|
| 642 | // Accumulate variance per if/pol/beam averaged over spectrum
|
---|
| 643 | // This method to get pos2 from itData.pos() is only valid
|
---|
| 644 | // because the spectral axis is the last one (so we can just
|
---|
| 645 | // copy the first nAXesSub positions out)
|
---|
| 646 |
|
---|
| 647 | pos = itData.pos().getFirst(nAxesSub);
|
---|
| 648 | sumSq(pos) += fac;
|
---|
| 649 | //
|
---|
| 650 | itData.next();
|
---|
| 651 | itMask.next();
|
---|
| 652 | }
|
---|
| 653 | } else if (wtType==TSYS) {
|
---|
| 654 | }
|
---|
| 655 |
|
---|
| 656 | // Accumulate sum of (possibly scaled) data
|
---|
| 657 |
|
---|
| 658 | sum += dataIn;
|
---|
| 659 |
|
---|
| 660 | // Accumulate Tsys, time, and interval
|
---|
| 661 |
|
---|
| 662 | tSysSum += tSys;
|
---|
| 663 | timeSum += time;
|
---|
| 664 | intSum += interval;
|
---|
| 665 | nAccum += 1;
|
---|
| 666 | }
|
---|
| 667 |
|
---|
[152] | 668 | SDMemTable* SDMath::localOperate (const SDMemTable& in, Float val, Bool doAll,
|
---|
| 669 | uInt what)
|
---|
| 670 | //
|
---|
| 671 | // what = 0 Multiply
|
---|
| 672 | // 1 Add
|
---|
[146] | 673 | {
|
---|
[152] | 674 | SDMemTable* pOut = new SDMemTable(in,False);
|
---|
| 675 | const Table& tOut = pOut->table();
|
---|
| 676 | ArrayColumn<Float> spec(tOut,"SPECTRA");
|
---|
[146] | 677 | //
|
---|
[152] | 678 | if (doAll) {
|
---|
| 679 | for (uInt i=0; i < tOut.nrow(); i++) {
|
---|
| 680 |
|
---|
| 681 | // Get
|
---|
| 682 |
|
---|
| 683 | MaskedArray<Float> marr(pOut->rowAsMaskedArray(i));
|
---|
| 684 |
|
---|
| 685 | // Operate
|
---|
| 686 |
|
---|
| 687 | if (what==0) {
|
---|
| 688 | marr *= val;
|
---|
| 689 | } else if (what==1) {
|
---|
| 690 | marr += val;
|
---|
| 691 | }
|
---|
| 692 |
|
---|
| 693 | // Put
|
---|
| 694 |
|
---|
| 695 | spec.put(i, marr.getArray());
|
---|
| 696 | }
|
---|
| 697 | } else {
|
---|
| 698 |
|
---|
| 699 | // Get cursor location
|
---|
| 700 |
|
---|
| 701 | IPosition start, end;
|
---|
| 702 | getCursorLocation (start, end, in);
|
---|
| 703 | //
|
---|
| 704 | for (uInt i=0; i < tOut.nrow(); i++) {
|
---|
| 705 |
|
---|
| 706 | // Get
|
---|
| 707 |
|
---|
| 708 | MaskedArray<Float> dataIn(pOut->rowAsMaskedArray(i));
|
---|
| 709 |
|
---|
| 710 | // Modify. More work than we would like to deal with the mask
|
---|
| 711 |
|
---|
| 712 | Array<Float>& values = dataIn.getRWArray();
|
---|
| 713 | Array<Bool> mask(dataIn.getMask());
|
---|
| 714 | //
|
---|
| 715 | Array<Float> values2 = values(start,end);
|
---|
| 716 | Array<Bool> mask2 = mask(start,end);
|
---|
| 717 | MaskedArray<Float> t(values2,mask2);
|
---|
| 718 | if (what==0) {
|
---|
| 719 | t *= val;
|
---|
| 720 | } else if (what==1) {
|
---|
| 721 | t += val;
|
---|
| 722 | }
|
---|
| 723 | values(start, end) = t.getArray(); // Write back into 'dataIn'
|
---|
| 724 |
|
---|
| 725 | // Put
|
---|
| 726 | spec.put(i, dataIn.getArray());
|
---|
| 727 | }
|
---|
| 728 | }
|
---|
| 729 | //
|
---|
[146] | 730 | return pOut;
|
---|
| 731 | }
|
---|
| 732 |
|
---|
| 733 |
|
---|
[152] | 734 |
|
---|
| 735 | void SDMath::getCursorLocation (IPosition& start, IPosition& end,
|
---|
| 736 | const SDMemTable& in)
|
---|
| 737 | {
|
---|
| 738 | const uInt nDim = 4;
|
---|
| 739 | const uInt i = in.getBeam();
|
---|
| 740 | const uInt j = in.getIF();
|
---|
| 741 | const uInt k = in.getPol();
|
---|
| 742 | const uInt n = in.nChan();
|
---|
| 743 | //
|
---|
[163] | 744 | start.resize(nDim);
|
---|
| 745 | start(0) = i;
|
---|
| 746 | start(1) = j;
|
---|
| 747 | start(2) = k;
|
---|
| 748 | start(3) = 0;
|
---|
[152] | 749 | //
|
---|
| 750 | end.resize(nDim);
|
---|
[163] | 751 | end(0) = i;
|
---|
| 752 | end(1) = j;
|
---|
| 753 | end(2) = k;
|
---|
| 754 | end(3) = n-1;
|
---|
[152] | 755 | }
|
---|
[163] | 756 |
|
---|
| 757 |
|
---|
| 758 | void SDMath::convertWeightString (WeightType& wtType, const std::string& weightStr)
|
---|
| 759 | {
|
---|
| 760 | String tStr(weightStr);
|
---|
| 761 | tStr.upcase();
|
---|
| 762 | if (tStr.contains(String("NONE"))) {
|
---|
| 763 | wtType = NONE;
|
---|
| 764 | } else if (tStr.contains(String("VAR"))) {
|
---|
| 765 | wtType = VAR;
|
---|
| 766 | } else if (tStr.contains(String("TSYS"))) {
|
---|
| 767 | wtType = TSYS;
|
---|
| 768 | throw (AipsError("T_sys weighting not yet implemented"));
|
---|
| 769 | } else {
|
---|
| 770 | throw (AipsError("Unrecognized weighting type"));
|
---|
| 771 | }
|
---|
| 772 | }
|
---|
| 773 |
|
---|
| 774 | void SDMath::putDataInSDC (SDContainer& sc, const Array<Float>& data,
|
---|
| 775 | const Array<Bool>& mask)
|
---|
| 776 | {
|
---|
| 777 | sc.putSpectrum(data);
|
---|
| 778 | //
|
---|
| 779 | Array<uChar> outflags(data.shape());
|
---|
| 780 | convertArray(outflags,!mask);
|
---|
| 781 | sc.putFlags(outflags);
|
---|
| 782 | }
|
---|
[165] | 783 |
|
---|
| 784 |
|
---|
| 785 | SDMemTable* SDMath::localAveragePol(const SDMemTable& in, const Vector<Bool>& mask)
|
---|
| 786 | //
|
---|
| 787 | // Average all polarizations together, weighted by variance
|
---|
| 788 | //
|
---|
| 789 | {
|
---|
| 790 | // WeightType wtType = NONE;
|
---|
| 791 | // convertWeightString (wtType, weight);
|
---|
| 792 |
|
---|
| 793 | const uInt nRows = in.nRow();
|
---|
| 794 | const uInt polAxis = 2; // Polarization axis
|
---|
| 795 | const uInt chanAxis = 3; // Spectrum axis
|
---|
| 796 |
|
---|
| 797 | // Create output Table and reshape number of polarizations
|
---|
| 798 |
|
---|
| 799 | Bool clear=True;
|
---|
| 800 | SDMemTable* pTabOut = new SDMemTable(in, clear);
|
---|
| 801 | SDHeader header = pTabOut->getSDHeader();
|
---|
| 802 | header.npol = 1;
|
---|
| 803 | pTabOut->putSDHeader(header);
|
---|
| 804 |
|
---|
| 805 | // Shape of input and output data
|
---|
| 806 |
|
---|
| 807 | const IPosition& shapeIn = in.rowAsMaskedArray(0u, False).shape();
|
---|
| 808 | IPosition shapeOut(shapeIn);
|
---|
| 809 | shapeOut(polAxis) = 1; // Average all polarizations
|
---|
| 810 | //
|
---|
| 811 | const uInt nChan = shapeIn(chanAxis);
|
---|
| 812 | const IPosition vecShapeOut(4,1,1,1,nChan); // A multi-dim form of a Vector shape
|
---|
| 813 | IPosition start(4), end(4);
|
---|
| 814 |
|
---|
| 815 | // Output arrays
|
---|
| 816 |
|
---|
| 817 | Array<Float> outData(shapeOut, 0.0);
|
---|
| 818 | Array<Bool> outMask(shapeOut, True);
|
---|
| 819 | const IPosition axes(2, 2, 3); // pol-channel plane
|
---|
| 820 | //
|
---|
| 821 | const Bool useMask = (mask.nelements() == shapeIn(chanAxis));
|
---|
| 822 |
|
---|
| 823 | // Loop over rows
|
---|
| 824 |
|
---|
| 825 | for (uInt iRow=0; iRow<nRows; iRow++) {
|
---|
| 826 |
|
---|
| 827 | // Get data for this row
|
---|
| 828 |
|
---|
| 829 | MaskedArray<Float> marr(in.rowAsMaskedArray(iRow));
|
---|
| 830 | Array<Float>& arr = marr.getRWArray();
|
---|
| 831 | const Array<Bool>& barr = marr.getMask();
|
---|
| 832 |
|
---|
| 833 | // Make iterators to iterate by pol-channel planes
|
---|
| 834 |
|
---|
| 835 | ReadOnlyArrayIterator<Float> itDataPlane(arr, axes);
|
---|
| 836 | ReadOnlyArrayIterator<Bool> itMaskPlane(barr, axes);
|
---|
| 837 |
|
---|
| 838 | // Accumulations
|
---|
| 839 |
|
---|
| 840 | Float fac = 1.0;
|
---|
| 841 | Vector<Float> vecSum(nChan,0.0);
|
---|
| 842 |
|
---|
| 843 | // Iterate through data by pol-channel planes
|
---|
| 844 |
|
---|
| 845 | while (!itDataPlane.pastEnd()) {
|
---|
| 846 |
|
---|
| 847 | // Iterate through plane by polarization and accumulate Vectors
|
---|
| 848 |
|
---|
| 849 | Vector<Float> t1(nChan); t1 = 0.0;
|
---|
| 850 | Vector<Bool> t2(nChan); t2 = True;
|
---|
| 851 | MaskedArray<Float> vecSum(t1,t2);
|
---|
| 852 | Float varSum = 0.0;
|
---|
| 853 | {
|
---|
| 854 | ReadOnlyVectorIterator<Float> itDataVec(itDataPlane.array(), 1);
|
---|
| 855 | ReadOnlyVectorIterator<Bool> itMaskVec(itMaskPlane.array(), 1);
|
---|
| 856 | while (!itDataVec.pastEnd()) {
|
---|
| 857 |
|
---|
| 858 | // Create MA of data & mask (optionally including OTF mask) and get variance
|
---|
| 859 |
|
---|
| 860 | if (useMask) {
|
---|
| 861 | const MaskedArray<Float> spec(itDataVec.vector(),mask&&itMaskVec.vector());
|
---|
| 862 | fac = 1.0 / variance(spec);
|
---|
| 863 | } else {
|
---|
| 864 | const MaskedArray<Float> spec(itDataVec.vector(),itMaskVec.vector());
|
---|
| 865 | fac = 1.0 / variance(spec);
|
---|
| 866 | }
|
---|
| 867 |
|
---|
| 868 | // Normalize spectrum (without OTF mask) and accumulate
|
---|
| 869 |
|
---|
| 870 | const MaskedArray<Float> spec(fac*itDataVec.vector(), itMaskVec.vector());
|
---|
| 871 | vecSum += spec;
|
---|
| 872 | varSum += fac;
|
---|
| 873 |
|
---|
| 874 | // Next
|
---|
| 875 |
|
---|
| 876 | itDataVec.next();
|
---|
| 877 | itMaskVec.next();
|
---|
| 878 | }
|
---|
| 879 | }
|
---|
| 880 |
|
---|
| 881 | // Normalize summed spectrum
|
---|
| 882 |
|
---|
| 883 | vecSum /= varSum;
|
---|
| 884 |
|
---|
| 885 | // FInd position in input data array. We are iterating by pol-channel
|
---|
| 886 | // plane so all that will change is beam and IF and that's what we want.
|
---|
| 887 |
|
---|
| 888 | IPosition pos = itDataPlane.pos();
|
---|
| 889 |
|
---|
| 890 | // Write out data. This is a bit messy. We have to reform the Vector
|
---|
| 891 | // accumulator into an Array of shape (1,1,1,nChan)
|
---|
| 892 |
|
---|
| 893 | start = pos;
|
---|
| 894 | end = pos;
|
---|
| 895 | end(chanAxis) = nChan-1;
|
---|
| 896 | outData(start,end) = vecSum.getArray().reform(vecShapeOut);
|
---|
| 897 | outMask(start,end) = vecSum.getMask().reform(vecShapeOut);
|
---|
| 898 |
|
---|
| 899 | // Step to next beam/IF combination
|
---|
| 900 |
|
---|
| 901 | itDataPlane.next();
|
---|
| 902 | itMaskPlane.next();
|
---|
| 903 | }
|
---|
| 904 |
|
---|
| 905 | // Generate output container and write it to output table
|
---|
| 906 |
|
---|
| 907 | SDContainer sc = in.getSDContainer();
|
---|
| 908 | sc.resize(shapeOut);
|
---|
| 909 | //
|
---|
| 910 | putDataInSDC (sc, outData, outMask);
|
---|
| 911 | pTabOut->putSDContainer(sc);
|
---|
| 912 | }
|
---|
| 913 | //
|
---|
| 914 | return pTabOut;
|
---|
| 915 | }
|
---|
[167] | 916 |
|
---|
| 917 | SDMemTable* SDMath::localBin (const SDMemTable& in, Int width)
|
---|
| 918 | {
|
---|
| 919 | SDHeader sh = in.getSDHeader();
|
---|
| 920 | SDMemTable* pTabOut = new SDMemTable(in, True);
|
---|
| 921 |
|
---|
| 922 | // Bin up SpectralCoordinates
|
---|
| 923 |
|
---|
| 924 | IPosition factors(1);
|
---|
| 925 | factors(0) = width;
|
---|
| 926 | for (uInt j=0; j<in.nCoordinates(); ++j) {
|
---|
| 927 | CoordinateSystem cSys;
|
---|
| 928 | cSys.addCoordinate(in.getCoordinate(j));
|
---|
| 929 | CoordinateSystem cSysBin =
|
---|
| 930 | CoordinateUtil::makeBinnedCoordinateSystem (factors, cSys, False);
|
---|
| 931 | //
|
---|
| 932 | SpectralCoordinate sCBin = cSysBin.spectralCoordinate(0);
|
---|
| 933 | pTabOut->setCoordinate(sCBin, j);
|
---|
| 934 | }
|
---|
| 935 |
|
---|
| 936 | // Use RebinLattice to find shape
|
---|
| 937 |
|
---|
| 938 | IPosition shapeIn(1,sh.nchan);
|
---|
| 939 | IPosition shapeOut = RebinLattice<Float>::rebinShape (shapeIn, factors);
|
---|
| 940 | sh.nchan = shapeOut(0);
|
---|
| 941 | pTabOut->putSDHeader(sh);
|
---|
| 942 |
|
---|
| 943 |
|
---|
| 944 | // Loop over rows and bin along channel axis
|
---|
| 945 |
|
---|
| 946 | const uInt axis = 3;
|
---|
| 947 | for (uInt i=0; i < in.nRow(); ++i) {
|
---|
| 948 | SDContainer sc = in.getSDContainer(i);
|
---|
| 949 | //
|
---|
| 950 | Array<Float> tSys(sc.getTsys()); // Get it out before sc changes shape
|
---|
| 951 |
|
---|
| 952 | // Bin up spectrum
|
---|
| 953 |
|
---|
| 954 | MaskedArray<Float> marr(in.rowAsMaskedArray(i));
|
---|
| 955 | MaskedArray<Float> marrout;
|
---|
| 956 | LatticeUtilities::bin(marrout, marr, axis, width);
|
---|
| 957 |
|
---|
| 958 | // Put back the binned data and flags
|
---|
| 959 |
|
---|
| 960 | IPosition ip2 = marrout.shape();
|
---|
| 961 | sc.resize(ip2);
|
---|
| 962 | //
|
---|
| 963 | putDataInSDC (sc, marrout.getArray(), marrout.getMask());
|
---|
| 964 |
|
---|
| 965 | // Bin up Tsys.
|
---|
| 966 |
|
---|
| 967 | Array<Bool> allGood(tSys.shape(),True);
|
---|
| 968 | MaskedArray<Float> tSysIn(tSys, allGood, True);
|
---|
| 969 | //
|
---|
| 970 | MaskedArray<Float> tSysOut;
|
---|
| 971 | LatticeUtilities::bin(tSysOut, tSysIn, axis, width);
|
---|
| 972 | sc.putTsys(tSysOut.getArray());
|
---|
| 973 | //
|
---|
| 974 | pTabOut->putSDContainer(sc);
|
---|
| 975 | }
|
---|
| 976 | return pTabOut;
|
---|
| 977 | }
|
---|
| 978 |
|
---|
| 979 |
|
---|