Line data Source code
1 : //# RegriddingTVI.h: This file contains the implementation of the RegriddingTVI class.
2 : //#
3 : //# CASA - Common Astronomy Software Applications (http://casa.nrao.edu/)
4 : //# Copyright (C) Associated Universities, Inc. Washington DC, USA 2011, All rights reserved.
5 : //# Copyright (C) European Southern Observatory, 2011, All rights reserved.
6 : //#
7 : //# This library is free software; you can redistribute it and/or
8 : //# modify it under the terms of the GNU Lesser General Public
9 : //# License as published by the Free software Foundation; either
10 : //# version 2.1 of the License, or (at your option) any later version.
11 : //#
12 : //# This library is distributed in the hope that it will be useful,
13 : //# but WITHOUT ANY WARRANTY, without even the implied warranty of
14 : //# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15 : //# Lesser General Public License for more details.
16 : //#
17 : //# You should have received a copy of the GNU Lesser General Public
18 : //# License along with this library; if not, write to the Free Software
19 : //# Foundation, Inc., 59 Temple Place, Suite 330, Boston,
20 : //# MA 02111-1307 USA
21 : //# $Id: $
22 :
23 : #include <mstransform/TVI/RegriddingTVI.h>
24 :
25 : using namespace casacore;
26 : namespace casa { //# NAMESPACE CASA - BEGIN
27 :
28 : namespace vi { //# NAMESPACE VI - BEGIN
29 :
30 : //////////////////////////////////////////////////////////////////////////
31 : // RegriddingTVI class
32 : //////////////////////////////////////////////////////////////////////////
33 :
34 : // -----------------------------------------------------------------------
35 : //
36 : // -----------------------------------------------------------------------
37 0 : RegriddingTVI::RegriddingTVI( ViImplementation2 * inputVii,
38 0 : const Record &configuration):
39 0 : FreqAxisTVI (inputVii)
40 : {
41 : // Frequency specification parameters
42 0 : nChan_p = -1;
43 0 : mode_p = String("channel"); // Options are: channel, frequency, velocity
44 0 : start_p = String("0");
45 0 : width_p = String("1"); // -1 means use all the input channels
46 0 : velocityType_p = String("radio"); // When mode is velocity options are: optical, radio
47 0 : restFrequency_p = String("");
48 0 : interpolationMethodPar_p = String("linear"); // Options are: nearest, linear, cubic, spline, fftshift
49 0 : outputReferenceFramePar_p = String(""); // Options are: LSRK, LSRD, BARY, GALACTO, LGROUP, CMB, GEO, or TOPO
50 0 : phaseCenterPar_p = new casac::variant("");
51 0 : regriddingMethod_p = linear;
52 :
53 : // Sub-cases
54 0 : refFrameTransformation_p = false;
55 0 : radialVelocityCorrection_p = false;
56 0 : fftShift_p = 0;
57 :
58 : // SPW-indexed maps
59 0 : weightFactorMap_p.clear();
60 0 : sigmaFactorMap_p.clear();
61 0 : inputOutputSpwMap_p.clear();
62 0 : freqTransEngineRowId_p = UINT_MAX;
63 :
64 : // Parse and check configuration parameters
65 : // Note: if a constructor finishes by throwing an exception, the memory
66 : // associated with the object itself is cleaned up — there is no memory leak.
67 0 : if (not parseConfiguration(configuration))
68 : {
69 0 : throw AipsError("Error parsing RegriddingTVI configuration");
70 : }
71 :
72 0 : initialize();
73 :
74 0 : return;
75 : }
76 :
77 : // -----------------------------------------------------------------------
78 : //
79 : // -----------------------------------------------------------------------
80 0 : Bool RegriddingTVI::parseConfiguration(const Record &configuration)
81 : {
82 0 : int exists = -1;
83 0 : Bool ret = true;
84 :
85 0 : exists = configuration.fieldNumber ("phasecenter");
86 0 : if (exists >= 0)
87 : {
88 : //If phase center is a simple numeric value then it is taken
89 : // as a FIELD_ID otherwise it is converted to a MDirection
90 0 : if( configuration.type(exists) == TpInt )
91 : {
92 0 : int fieldIdForPhaseCenter = -1;
93 0 : configuration.get (exists, fieldIdForPhaseCenter);
94 0 : logger_p << LogIO::NORMAL << LogOrigin("RegriddingTVI", __FUNCTION__)
95 0 : << "Field Id for phase center is " << fieldIdForPhaseCenter << LogIO::POST;
96 0 : if (phaseCenterPar_p) delete phaseCenterPar_p;
97 0 : phaseCenterPar_p = new casac::variant((long)fieldIdForPhaseCenter);
98 : }
99 : else
100 : {
101 0 : String phaseCenter("");
102 0 : configuration.get (exists, phaseCenter);
103 0 : logger_p << LogIO::NORMAL << LogOrigin("RegriddingTVI", __FUNCTION__)
104 0 : << "Phase center is " << phaseCenter << LogIO::POST;
105 0 : if (phaseCenterPar_p) delete phaseCenterPar_p;
106 0 : phaseCenterPar_p = new casac::variant(phaseCenter);
107 : }
108 : }
109 :
110 0 : exists = -1;
111 0 : exists = configuration.fieldNumber ("restfreq");
112 0 : if (exists >= 0)
113 : {
114 0 : configuration.get (exists, restFrequency_p);
115 0 : if (!restFrequency_p.empty())
116 : {
117 0 : logger_p << LogIO::NORMAL << LogOrigin("RegriddingTVI", __FUNCTION__)
118 0 : << "Rest frequency is " << restFrequency_p << LogIO::POST;
119 : }
120 : }
121 :
122 0 : exists = configuration.fieldNumber ("outframe");
123 0 : if (exists >= 0)
124 : {
125 0 : configuration.get (exists, outputReferenceFramePar_p);
126 0 : logger_p << LogIO::NORMAL << LogOrigin("RegriddingTVI", __FUNCTION__)
127 0 : << "Output reference frame is " << outputReferenceFramePar_p << LogIO::POST;
128 : }
129 :
130 0 : exists = -1;
131 0 : exists = configuration.fieldNumber ("interpolation");
132 0 : if (exists >= 0)
133 : {
134 0 : configuration.get (exists, interpolationMethodPar_p);
135 :
136 0 : if (interpolationMethodPar_p.contains("nearest"))
137 : {
138 0 : regriddingMethod_p = nearestNeighbour;
139 : }
140 0 : else if (interpolationMethodPar_p.contains("linear"))
141 : {
142 0 : regriddingMethod_p = linear;
143 : }
144 0 : else if (interpolationMethodPar_p.contains("cubic"))
145 : {
146 0 : regriddingMethod_p = cubic;
147 : }
148 0 : else if (interpolationMethodPar_p.contains("spline"))
149 : {
150 0 : regriddingMethod_p = spline;
151 : }
152 0 : else if (interpolationMethodPar_p.contains("fftshift"))
153 : {
154 0 : regriddingMethod_p = fftshift;
155 : }
156 : else
157 : {
158 0 : logger_p << LogIO::SEVERE << LogOrigin("RegriddingTVI", __FUNCTION__)
159 0 : << "Interpolation method " << interpolationMethodPar_p << " not available " << LogIO::POST;
160 : }
161 :
162 0 : logger_p << LogIO::NORMAL << LogOrigin("RegriddingTVI", __FUNCTION__)
163 0 : << "Interpolation method is " << interpolationMethodPar_p << LogIO::POST;
164 : }
165 : else
166 : {
167 0 : regriddingMethod_p = linear;
168 : }
169 :
170 0 : exists = -1;
171 0 : exists = configuration.fieldNumber ("mode");
172 0 : if (exists >= 0)
173 : {
174 0 : configuration.get (exists, mode_p);
175 0 : logger_p << LogIO::NORMAL << LogOrigin("RegriddingTVI", __FUNCTION__)
176 0 : << "Mode is " << mode_p<< LogIO::POST;
177 :
178 0 : if ((mode_p == "frequency") or (mode_p == "velocity"))
179 : {
180 0 : start_p = String("");
181 0 : width_p = String("");
182 : }
183 : }
184 :
185 0 : exists = -1;
186 0 : exists = configuration.fieldNumber ("nchan");
187 0 : if (exists >= 0)
188 : {
189 0 : configuration.get (exists, nChan_p);
190 0 : logger_p << LogIO::NORMAL << LogOrigin("RegriddingTVI", __FUNCTION__)
191 0 : << "Number of channels is " << nChan_p<< LogIO::POST;
192 : }
193 :
194 0 : exists = -1;
195 0 : exists = configuration.fieldNumber ("start");
196 0 : if (exists >= 0)
197 : {
198 0 : configuration.get (exists, start_p);
199 0 : logger_p << LogIO::NORMAL << LogOrigin("RegriddingTVI", __FUNCTION__)
200 0 : << "Start is " << start_p << LogIO::POST;
201 : }
202 :
203 0 : exists = -1;
204 0 : exists = configuration.fieldNumber ("width");
205 0 : if (exists >= 0)
206 : {
207 0 : configuration.get (exists, width_p);
208 0 : logger_p << LogIO::NORMAL << LogOrigin("RegriddingTVI", __FUNCTION__)
209 0 : << "Width is " << width_p << LogIO::POST;
210 : }
211 :
212 0 : exists = -1;
213 0 : exists = configuration.fieldNumber ("veltype");
214 0 : if ((exists >= 0) and (mode_p == "velocity"))
215 : {
216 0 : configuration.get (exists, velocityType_p);
217 0 : logger_p << LogIO::NORMAL << LogOrigin("RegriddingTVI", __FUNCTION__)
218 0 : << "Velocity type is " << velocityType_p << LogIO::POST;
219 : }
220 :
221 0 : return ret;
222 : }
223 :
224 : // -----------------------------------------------------------------------
225 : //
226 : // -----------------------------------------------------------------------
227 0 : void RegriddingTVI::initialize()
228 : {
229 : // Determine input reference frame from the first row in the SPW sub-table of the output (selected) MS
230 0 : MSSpectralWindow spwTable = getVii()->ms().spectralWindow();
231 0 : MSSpWindowColumns spwCols(spwTable);
232 0 : inputReferenceFrame_p = MFrequency::castType(spwCols.measFreqRef()(0));
233 :
234 : // Parse output reference frame
235 0 : refFrameTransformation_p = true;
236 0 : radialVelocityCorrection_p = false;
237 0 : if(outputReferenceFramePar_p.empty())
238 : {
239 0 : outputReferenceFrame_p = inputReferenceFrame_p;
240 : }
241 : // CAS-6778: Support for new ref. frame SOURCE that requires radial velocity correction
242 0 : else if (outputReferenceFramePar_p == "SOURCE")
243 : {
244 0 : outputReferenceFrame_p = MFrequency::GEO;
245 0 : radialVelocityCorrection_p = true;
246 : }
247 0 : else if(!MFrequency::getType(outputReferenceFrame_p, outputReferenceFramePar_p))
248 : {
249 0 : logger_p << LogIO::SEVERE << LogOrigin("RegriddingTVI", __FUNCTION__)
250 0 : << "Problem parsing output reference frame:" << outputReferenceFramePar_p << LogIO::EXCEPTION;
251 : }
252 :
253 0 : if (outputReferenceFrame_p == inputReferenceFrame_p) {
254 0 : refFrameTransformation_p = false;
255 : }
256 :
257 :
258 : // Determine observatory position from the first row in the observation sub-table of the output (selected) MS
259 0 : MSObservation observationTable = getVii()->ms().observation();
260 0 : MSObservationColumns observationCols(observationTable);
261 0 : String observatoryName = observationCols.telescopeName()(0);
262 0 : MeasTable::Observatory(observatoryPosition_p,observatoryName);
263 :
264 : // jagonzal: This conversion is needed only for cosmetic reasons
265 : // observatoryPosition_p=MPosition::Convert(observatoryPosition_p, MPosition::ITRF)();
266 :
267 : // Determine observation time from the first row in the selected MS
268 0 : selectedInputMsCols_p = new MSColumns(getVii()->ms());
269 0 : referenceTime_p = selectedInputMsCols_p->timeMeas()(0);
270 :
271 : // Access FIELD cols to get phase center and radial velocity
272 0 : MSField field = getVii()->ms().field();
273 0 : inputMSFieldCols_p = new MSFieldColumns(field);
274 :
275 : // Determine phase center
276 0 : if (phaseCenterPar_p->type() == casac::variant::INT)
277 : {
278 0 : Int fieldIdForPhaseCenter = phaseCenterPar_p->toInt();
279 :
280 0 : if (fieldIdForPhaseCenter >= (Int)inputMSFieldCols_p->nrow() || fieldIdForPhaseCenter < 0)
281 : {
282 0 : logger_p << LogIO::SEVERE << LogOrigin("RegriddingTVI", __FUNCTION__)
283 0 : << "Selected FIELD_ID to determine phase center does not exist " << LogIO::POST;
284 : }
285 : else
286 : {
287 : // CAS-6778: Support for new ref. frame SOURCE that requires radial velocity correction
288 0 : if (radialVelocityCorrection_p)
289 : {
290 0 : radialVelocity_p = inputMSFieldCols_p->radVelMeas(fieldIdForPhaseCenter, referenceTime_p.get("s").getValue());
291 0 : phaseCenter_p = inputMSFieldCols_p->phaseDirMeas(fieldIdForPhaseCenter,referenceTime_p.get("s").getValue());
292 : }
293 : else
294 : {
295 0 : phaseCenter_p = inputMSFieldCols_p->phaseDirMeasCol()(fieldIdForPhaseCenter)(IPosition(1,0));
296 : }
297 : }
298 : }
299 : else
300 : {
301 0 : String phaseCenter = phaseCenterPar_p->toString(true);
302 :
303 : // Determine phase center from the first row in the FIELD sub-table of the output (selected) MS
304 0 : if (phaseCenter.empty())
305 : {
306 : // CAS-6778: Support for new ref. frame SOURCE that requires radial velocity correction
307 0 : if (radialVelocityCorrection_p)
308 : {
309 0 : radialVelocity_p = inputMSFieldCols_p->radVelMeas(0, referenceTime_p.get("s").getValue());
310 0 : phaseCenter_p = inputMSFieldCols_p->phaseDirMeas(0,referenceTime_p.get("s").getValue());
311 : }
312 : else
313 : {
314 0 : phaseCenter_p = inputMSFieldCols_p->phaseDirMeasCol()(0)(IPosition(1,0));
315 : }
316 : }
317 : // Parse phase center
318 : else
319 : {
320 0 : if(!casaMDirection(phaseCenter, phaseCenter_p))
321 : {
322 0 : logger_p << LogIO::SEVERE << LogOrigin("MSTransformManager", __FUNCTION__)
323 0 : << "Cannot interpret phase center " << phaseCenter << LogIO::POST;
324 0 : return;
325 : }
326 : }
327 : }
328 :
329 0 : if (radialVelocityCorrection_p && (radialVelocity_p.getRef().getType() != MRadialVelocity::GEO))
330 : {
331 0 : logger_p << LogIO::SEVERE << LogOrigin("RegriddingTVI", __FUNCTION__)
332 : << "Cannot perform radial velocity correction with ephemerides of type "
333 0 : << MRadialVelocity::showType(radialVelocity_p.getRef().getType()) << ".\nType needs to be GEO."
334 0 : << LogIO::EXCEPTION;
335 : }
336 :
337 0 : return;
338 : }
339 :
340 : // -----------------------------------------------------------------------
341 : //
342 : // -----------------------------------------------------------------------
343 0 : void RegriddingTVI::origin()
344 : {
345 : // Drive underlying ViImplementation2
346 0 : getVii()->origin();
347 :
348 : // Define output grid for this chunk (also defines shape)
349 0 : initFrequencyGrid();
350 :
351 : // Define the shapes in the VB2
352 0 : configureShapes();
353 :
354 : // Synchronize own VisBuffer
355 0 : configureNewSubchunk();
356 :
357 0 : return;
358 : }
359 :
360 : // -----------------------------------------------------------------------
361 : //
362 : // -----------------------------------------------------------------------
363 0 : Int RegriddingTVI::getReportingFrameOfReference () const
364 : {
365 0 : return outputReferenceFrame_p;
366 : }
367 :
368 : // -----------------------------------------------------------------------
369 : //
370 : // -----------------------------------------------------------------------
371 0 : void RegriddingTVI::initFrequencyGrid()
372 : {
373 : // Get input VisBuffer
374 0 : VisBuffer2 *vb = getVii()->getVisBuffer();
375 :
376 : // Check if frequency grid has to be initialized for this SPW
377 0 : Int spwId = vb->spectralWindows()(0);
378 0 : if (inputOutputSpwMap_p.find(spwId) == inputOutputSpwMap_p.end())
379 : {
380 : // Get input frequencies in reporting frame of reference
381 : // i.e. as they are stored in the SPW sub-table
382 0 : Vector<Double> inputFreq = vb->getFrequencies(0);
383 0 : Vector<Double> inputWidth(inputFreq.size(),inputFreq(1)-inputFreq(0));
384 :
385 : // Declare output variables
386 : Double weightScale;
387 0 : Vector<Double> outputFreq;
388 0 : Vector<Double> outputWidth;
389 :
390 : // Use calcChanFreqs to change reference frame and regrid
391 0 : MFrequency::Types inputRefFrame = static_cast<MFrequency::Types> (getVii()->getReportingFrameOfReference());
392 0 : Bool ret = MSTransformRegridder::calcChanFreqs( logger_p,
393 : outputFreq,
394 : outputWidth,
395 : weightScale,
396 : inputFreq,
397 : inputWidth,
398 0 : phaseCenter_p,
399 : inputRefFrame,
400 0 : referenceTime_p,
401 0 : observatoryPosition_p,
402 0 : mode_p,
403 : nChan_p,
404 0 : start_p,
405 0 : width_p,
406 0 : restFrequency_p,
407 0 : outputReferenceFramePar_p,
408 0 : velocityType_p,
409 : true, // verbose
410 0 : radialVelocity_p
411 : );
412 :
413 0 : if (not ret)
414 : {
415 0 : throw AipsError("Error calculating output grid");
416 : }
417 :
418 : // Add input-output SPW pair to map
419 0 : spwInfo inputSpw(inputFreq,inputWidth);
420 0 : spwInfo outputSpw(outputFreq,outputWidth);
421 0 : inputOutputSpwMap_p[spwId] = std::make_pair(inputSpw,outputSpw);
422 :
423 : // Store weight/sigma factors
424 0 : weightFactorMap_p[spwId] = weightScale;
425 0 : sigmaFactorMap_p[spwId] = 1 / sqrt(weightScale);
426 :
427 : // Populate nchan input-output maps
428 0 : spwOutChanNumMap_p[spwId] = outputSpw.NUM_CHAN;
429 : }
430 :
431 0 : return;
432 : }
433 :
434 : // -----------------------------------------------------------------------
435 : //
436 : // -----------------------------------------------------------------------
437 0 : void RegriddingTVI::initFrequencyTransformationEngine() const
438 : {
439 : // Get input VisBuffer
440 0 : VisBuffer2 *vb = getVii()->getVisBuffer();
441 :
442 : // Check if frequency transformation engine has to be re-constructed
443 0 : auto rowId = vb->rowIds()[0];
444 0 : if (freqTransEngineRowId_p != rowId)
445 : {
446 : // Mark this rowId as the current one
447 0 : freqTransEngineRowId_p = rowId;
448 :
449 : // Get spwId to access input-output SPW map
450 0 : Int spwId = vb->spectralWindows()[0];
451 :
452 : // Get current time
453 0 : ScalarMeasColumn<MEpoch> mainTimeMeasCol = selectedInputMsCols_p->timeMeas();
454 0 : MEpoch currentRowTime = mainTimeMeasCol(rowId);
455 :
456 : // Check if radial velocity correction if necessary
457 0 : MDoppler radVelCorr;
458 0 : MDirection inputFieldDirection;
459 0 : Bool radVelSignificant = false;
460 0 : if (radialVelocityCorrection_p && inputMSFieldCols_p->needInterTime(vb->fieldId()(0)))
461 : {
462 0 : MRadialVelocity mRV = inputMSFieldCols_p->radVelMeas(vb->fieldId()(0),vb->time()(0));
463 0 : Quantity mrv = mRV.get("m/s");
464 0 : Quantity offsetMrv = radialVelocity_p.get("m/s"); // the radvel by which the out SPW def was shifted
465 0 : radVelCorr = MDoppler(mrv-(Quantity(2.)*offsetMrv));
466 0 : if (fabs(mrv.getValue()) > 1E-6) radVelSignificant = true;
467 :
468 0 : inputFieldDirection = inputMSFieldCols_p->phaseDirMeas(vb->fieldId()(0), vb->time()(0));
469 : }
470 : else
471 : {
472 0 : inputFieldDirection = vb->phaseCenter();
473 : }
474 :
475 : // Get input Ref. Frame (can be different for each SPWs)
476 0 : MFrequency::Types inputRefFrame = static_cast<MFrequency::Types> (getVii()->getReportingFrameOfReference());
477 :
478 : // Construct reference frame transformation engine
479 : MFrequency::Ref inputFrameRef = MFrequency::Ref(inputRefFrame,
480 0 : MeasFrame(inputFieldDirection,
481 : observatoryPosition_p,
482 0 : currentRowTime));
483 :
484 0 : MFrequency::Ref outputFrameRef = MFrequency::Ref(outputReferenceFrame_p,
485 0 : MeasFrame(phaseCenter_p,
486 : observatoryPosition_p,
487 0 : referenceTime_p));
488 :
489 : // Calculate new frequencies (also for FFT mode)
490 0 : freqTransEngine_p = MFrequency::Convert(Hz, inputFrameRef, outputFrameRef);
491 :
492 0 : for(uInt chan_idx=0; chan_idx<inputOutputSpwMap_p[spwId].first.CHAN_FREQ.size(); chan_idx++)
493 : {
494 0 : inputOutputSpwMap_p[spwId].first.CHAN_FREQ_aux[chan_idx] =
495 0 : freqTransEngine_p(inputOutputSpwMap_p[spwId].first.CHAN_FREQ[chan_idx]).get(Hz).getValue();
496 : }
497 :
498 : // Apply radial velocity correction if necessary
499 0 : if (radVelSignificant)
500 : {
501 0 : inputOutputSpwMap_p[spwId].first.CHAN_FREQ_aux =
502 0 : radVelCorr.shiftFrequency(inputOutputSpwMap_p[spwId].first.CHAN_FREQ_aux);
503 : }
504 :
505 : // Calculate FFT shift if necessary
506 0 : if (regriddingMethod_p == fftshift)
507 : {
508 0 : uInt centralChan = inputOutputSpwMap_p[spwId].first.CHAN_FREQ.size()/2;
509 0 : Double absoluteShift = inputOutputSpwMap_p[spwId].first.CHAN_FREQ_aux[centralChan]
510 0 : -inputOutputSpwMap_p[spwId].first.CHAN_FREQ[centralChan];
511 :
512 0 : Double chanWidth = inputOutputSpwMap_p[spwId].second.CHAN_FREQ[1]
513 0 : - inputOutputSpwMap_p[spwId].second.CHAN_FREQ[0];
514 :
515 0 : Double bandwidth = inputOutputSpwMap_p[spwId].second.CHAN_FREQ[inputOutputSpwMap_p[spwId].second.NUM_CHAN-1]
516 0 : - inputOutputSpwMap_p[spwId].second.CHAN_FREQ[0];
517 :
518 0 : bandwidth += chanWidth;
519 :
520 0 : fftShift_p = - absoluteShift / bandwidth;
521 : }
522 : }
523 :
524 0 : return;
525 : }
526 :
527 : // -----------------------------------------------------------------------
528 : //
529 : // -----------------------------------------------------------------------
530 0 : Vector<Double> RegriddingTVI::getFrequencies ( Double time,
531 : Int frameOfReference,
532 : Int spectralWindowId,
533 : Int msId) const
534 : {
535 0 : if (frameOfReference == outputReferenceFrame_p)
536 : {
537 0 : return inputOutputSpwMap_p[spectralWindowId].second.CHAN_FREQ;
538 : }
539 : else
540 : {
541 : // Get frequencies from input TVI
542 0 : Vector<Double> inputFreq = getVii()->getFrequencies(time,frameOfReference,spectralWindowId,msId);
543 :
544 0 : Vector<Double> inputWidth(inputFreq.size(),0);
545 0 : for (uInt chan_i=0;chan_i<inputFreq.size()-1;chan_i++)
546 : {
547 0 : inputWidth(chan_i) = inputFreq(chan_i+1)-inputFreq(chan_i);
548 : }
549 0 : inputWidth(inputFreq.size()-1) = inputWidth(inputFreq.size()-2);
550 :
551 : // Declare output variables
552 : Double weightScale;
553 0 : Vector<Double> outputFreq;
554 0 : Vector<Double> outputWidth;
555 :
556 : // Use calcChanFreqs for re-gridding only (output frame = inputFrame)
557 0 : MFrequency::Types frameOfReferenceEnum = static_cast<MFrequency::Types> (frameOfReference);
558 0 : String frameOfReferenceStr = MFrequency::showType(frameOfReferenceEnum);
559 0 : MSTransformRegridder::calcChanFreqs( logger_p,
560 : outputFreq,
561 : outputWidth,
562 : weightScale,
563 : inputFreq,
564 : inputWidth,
565 0 : phaseCenter_p,
566 : frameOfReferenceEnum,
567 0 : referenceTime_p,
568 0 : observatoryPosition_p,
569 0 : mode_p,
570 0 : nChan_p,
571 0 : start_p,
572 0 : width_p,
573 0 : restFrequency_p,
574 : frameOfReferenceStr,
575 0 : velocityType_p,
576 : false, // verbose
577 0 : radialVelocity_p
578 : );
579 :
580 0 : return outputFreq;
581 : }
582 :
583 : }
584 :
585 : // -----------------------------------------------------------------------
586 : //
587 : // -----------------------------------------------------------------------
588 0 : void RegriddingTVI::flag(Cube<Bool>& flagCube) const
589 : {
590 : // Get input VisBuffer and SPW
591 0 : VisBuffer2 *vb = getVii()->getVisBuffer();
592 0 : Int inputSPW = vb->spectralWindows()(0);
593 :
594 : // Configure Transformation Engine
595 0 : initFrequencyTransformationEngine();
596 :
597 : // Reshape output data before passing it to the DataCubeHolder
598 0 : flagCube.resize(getVisBuffer()->getShape(),false);
599 :
600 : // Gather input data
601 0 : DataCubeMap inputData;
602 0 : DataCubeHolder<Bool> inputFlagCubeHolder(vb->flagCube());
603 0 : inputData.add(MS::FLAG,inputFlagCubeHolder);
604 :
605 : // Gather output data
606 0 : DataCubeMap outputData;
607 0 : DataCubeHolder<Bool> outputFlagCubeHolder(flagCube);
608 0 : outputData.add(MS::FLAG,outputFlagCubeHolder);
609 :
610 : // Configure kernel
611 0 : if (regriddingMethod_p == fftshift)
612 : {
613 0 : DataFFTKernel<Float> kernel(regriddingMethod_p);
614 0 : RegriddingTransformEngine<Float> transformer(&kernel,&inputData,&outputData);
615 0 : transformFreqAxis2(vb->getShape(),transformer);
616 : }
617 : else
618 : {
619 0 : Vector<Double> *inputFreq = &(inputOutputSpwMap_p[inputSPW].first.CHAN_FREQ_aux);
620 0 : Vector<Double> *outputFreq = &(inputOutputSpwMap_p[inputSPW].second.CHAN_FREQ);
621 0 : DataInterpolationKernel<Float> kernel(regriddingMethod_p,inputFreq,outputFreq);
622 0 : RegriddingTransformEngine<Float> transformer(&kernel,&inputData,&outputData);
623 0 : transformFreqAxis2(vb->getShape(),transformer);
624 : }
625 :
626 0 : return;
627 : }
628 :
629 : // -----------------------------------------------------------------------
630 : //
631 : // -----------------------------------------------------------------------
632 0 : void RegriddingTVI::floatData (Cube<Float> & vis) const
633 : {
634 0 : transformDataCube(getVii()->getVisBuffer()->visCubeFloat(),vis);
635 :
636 0 : return;
637 : }
638 :
639 : // -----------------------------------------------------------------------
640 : //
641 : // -----------------------------------------------------------------------
642 0 : void RegriddingTVI::visibilityObserved (Cube<Complex> & vis) const
643 : {
644 0 : transformDataCube(getVii()->getVisBuffer()->visCube(),vis);
645 :
646 0 : return;
647 : }
648 :
649 : // -----------------------------------------------------------------------
650 : //
651 : // -----------------------------------------------------------------------
652 0 : void RegriddingTVI::visibilityCorrected (Cube<Complex> & vis) const
653 : {
654 0 : transformDataCube(getVii()->getVisBuffer()->visCubeCorrected(),vis);
655 :
656 0 : return;
657 : }
658 :
659 : // -----------------------------------------------------------------------
660 : //
661 : // -----------------------------------------------------------------------
662 0 : void RegriddingTVI::visibilityModel (Cube<Complex> & vis) const
663 : {
664 0 : transformDataCube(getVii()->getVisBuffer()->visCubeModel(),vis);
665 :
666 0 : return;
667 : }
668 :
669 : // -----------------------------------------------------------------------
670 : //
671 : // -----------------------------------------------------------------------
672 0 : void RegriddingTVI::weightSpectrum(Cube<Float> &weightSp) const
673 : {
674 : // Get input VisBuffer and SPW
675 0 : VisBuffer2 *vb = getVii()->getVisBuffer();
676 0 : Int inputSPW = vb->spectralWindows()(0);
677 :
678 : // Transform data
679 0 : transformDataCube(getVii()->getVisBuffer()->weightSpectrum(),weightSp);
680 :
681 : // Apply scaling factor on weights
682 0 : weightSp *= weightFactorMap_p[inputSPW];
683 :
684 :
685 0 : return;
686 : }
687 :
688 : // -----------------------------------------------------------------------
689 : //
690 : // -----------------------------------------------------------------------
691 0 : void RegriddingTVI::sigmaSpectrum (Cube<Float> &sigmaSp) const
692 : {
693 : // Get input VisBuffer and SPW
694 0 : VisBuffer2 *vb = getVii()->getVisBuffer();
695 0 : Int inputSPW = vb->spectralWindows()(0);
696 :
697 : // Transform data
698 0 : transformDataCube(getVii()->getVisBuffer()->sigmaSpectrum(),sigmaSp);
699 :
700 : // Apply scaling factor on weights
701 0 : sigmaSp *= sigmaFactorMap_p[inputSPW];
702 :
703 0 : return;
704 : }
705 :
706 : // -----------------------------------------------------------------------
707 : //
708 : // -----------------------------------------------------------------------
709 0 : template<class T> void RegriddingTVI::transformDataCube( const Cube<T> &inputVis,
710 : Cube<T> &outputVis) const
711 : {
712 : // Get input VisBuffer and SPW
713 0 : VisBuffer2 *vb = getVii()->getVisBuffer();
714 0 : Int inputSPW = vb->spectralWindows()(0);
715 :
716 : // Configure Transformation Engine
717 0 : initFrequencyTransformationEngine();
718 :
719 : // Reshape output data before passing it to the DataCubeHolder
720 0 : outputVis.resize(getVisBuffer()->getShape(),false);
721 :
722 : // Gather input data
723 0 : DataCubeMap inputData;
724 0 : DataCubeHolder<Bool> inputFlagCubeHolder(vb->flagCube());
725 0 : DataCubeHolder<T> inputVisCubeHolder(inputVis);
726 0 : inputData.add(MS::FLAG,inputFlagCubeHolder);
727 0 : inputData.add(MS::DATA,inputVisCubeHolder);
728 :
729 : // Gather output data
730 0 : DataCubeMap outputData;
731 0 : DataCubeHolder<T> outputVisCubeHolder(outputVis);
732 0 : outputData.add(MS::DATA,outputVisCubeHolder);
733 :
734 : // Configure kernel
735 0 : if (regriddingMethod_p == fftshift)
736 : {
737 0 : DataFFTKernel<T> kernel(regriddingMethod_p);
738 0 : RegriddingTransformEngine<T> transformer(&kernel,&inputData,&outputData);
739 0 : transformFreqAxis2(vb->getShape(),transformer);
740 : }
741 : else
742 : {
743 0 : Vector<Double> *inputFreq = &(inputOutputSpwMap_p[inputSPW].first.CHAN_FREQ_aux);
744 0 : Vector<Double> *outputFreq = &(inputOutputSpwMap_p[inputSPW].second.CHAN_FREQ);
745 0 : DataInterpolationKernel<T> kernel(regriddingMethod_p,inputFreq,outputFreq);
746 0 : RegriddingTransformEngine<T> transformer(&kernel,&inputData,&outputData);
747 0 : transformFreqAxis2(vb->getShape(),transformer);
748 : }
749 :
750 0 : return;
751 : }
752 :
753 : //////////////////////////////////////////////////////////////////////////
754 : // RegriddingTVIFactory class
755 : //////////////////////////////////////////////////////////////////////////
756 :
757 : // -----------------------------------------------------------------------
758 : //
759 : // -----------------------------------------------------------------------
760 0 : RegriddingTVIFactory::RegriddingTVIFactory (Record &configuration,
761 0 : ViImplementation2 *inputVii)
762 : {
763 0 : inputVii_p = inputVii;
764 0 : configuration_p = configuration;
765 0 : }
766 :
767 : // -----------------------------------------------------------------------
768 : //
769 : // -----------------------------------------------------------------------
770 0 : vi::ViImplementation2 * RegriddingTVIFactory::createVi(VisibilityIterator2 *) const
771 : {
772 0 : return new RegriddingTVI(inputVii_p,configuration_p);
773 : }
774 :
775 : // -----------------------------------------------------------------------
776 : //
777 : // -----------------------------------------------------------------------
778 :
779 0 : vi::ViImplementation2 * RegriddingTVIFactory::createVi() const
780 : {
781 0 : return new RegriddingTVI(inputVii_p,configuration_p);
782 : }
783 :
784 : //////////////////////////////////////////////////////////////////////////
785 : // RegriddingTVILayerFactory class
786 : //////////////////////////////////////////////////////////////////////////
787 :
788 0 : RegriddingTVILayerFactory::RegriddingTVILayerFactory(Record &configuration) :
789 : ViiLayerFactory(),
790 0 : configuration_p(configuration)
791 0 : {}
792 :
793 : ViImplementation2*
794 0 : RegriddingTVILayerFactory::createInstance(ViImplementation2* vii0) const
795 : {
796 : // Make the RegriddingTVi2, using supplied ViImplementation2, and return it
797 0 : ViImplementation2 *vii = new RegriddingTVI(vii0,configuration_p);
798 0 : return vii;
799 : }
800 :
801 : //////////////////////////////////////////////////////////////////////////
802 : // RegriddingTransformEngine class
803 : //////////////////////////////////////////////////////////////////////////
804 :
805 : // -----------------------------------------------------------------------
806 : //
807 : // -----------------------------------------------------------------------
808 0 : template<class T> RegriddingTransformEngine<T>::RegriddingTransformEngine( RegriddingKernel<T> *kernel,
809 : DataCubeMap *inputData,
810 : DataCubeMap *outputData):
811 : FreqAxisTransformEngine2<T>(inputData,
812 0 : outputData)
813 : {
814 :
815 0 : regriddingKernel_p = kernel;
816 0 : }
817 :
818 : // -----------------------------------------------------------------------
819 : //
820 : // -----------------------------------------------------------------------
821 0 : template<class T> void RegriddingTransformEngine<T>::transform()
822 : {
823 0 : regriddingKernel_p->kernel(inputData_p,outputData_p);
824 :
825 0 : return;
826 : }
827 :
828 : //////////////////////////////////////////////////////////////////////////
829 : // RegriddingKernel class
830 : //////////////////////////////////////////////////////////////////////////
831 :
832 : // -----------------------------------------------------------------------
833 : //
834 : // -----------------------------------------------------------------------
835 0 : template<class T> RegriddingKernel<T>::RegriddingKernel()
836 : {
837 0 : inputDummyFlagVectorInitialized_p = false;
838 0 : outputDummyFlagVectorInitialized_p = false;
839 0 : inputDummyDataVectorInitialized_p = false;
840 0 : outputDummyDataVectorInitialized_p = false;
841 0 : }
842 :
843 : // -----------------------------------------------------------------------
844 : //
845 : // -----------------------------------------------------------------------
846 0 : template<class T> Vector<Bool>& RegriddingKernel<T>::getInputFlagVector(DataCubeMap *inputData)
847 : {
848 0 : if (inputData->present(MS::FLAG))
849 : {
850 0 : return inputData->getVector<Bool>(MS::FLAG);
851 : }
852 0 : else if (not inputDummyFlagVectorInitialized_p)
853 : {
854 0 : inputDummyFlagVectorInitialized_p = true;
855 0 : inputDummyFlagVector_p.resize(inputData->getVectorShape()(0),false);
856 : }
857 :
858 0 : return inputDummyFlagVector_p;
859 : }
860 :
861 : // -----------------------------------------------------------------------
862 : //
863 : // -----------------------------------------------------------------------
864 0 : template<class T> Vector<Bool>& RegriddingKernel<T>::getOutputFlagVector(DataCubeMap *outputData)
865 : {
866 0 : if (outputData->present(MS::FLAG))
867 : {
868 0 : return outputData->getVector<Bool>(MS::FLAG);
869 : }
870 0 : else if (not outputDummyFlagVectorInitialized_p)
871 : {
872 0 : outputDummyFlagVectorInitialized_p = true;
873 0 : outputDummyFlagVector_p.resize(outputData->getVectorShape()(0),false);
874 : }
875 :
876 0 : return outputDummyFlagVector_p;
877 : }
878 :
879 : // -----------------------------------------------------------------------
880 : //
881 : // -----------------------------------------------------------------------
882 0 : template<class T> Vector<T>& RegriddingKernel<T>::getInputDataVector(DataCubeMap *inputData)
883 : {
884 0 : if (inputData->present(MS::DATA))
885 : {
886 0 : return inputData->getVector<T>(MS::DATA);
887 : }
888 0 : else if (not inputDummyDataVectorInitialized_p)
889 : {
890 0 : inputDummyDataVectorInitialized_p = true;
891 0 : inputDummyDataVector_p.resize(inputData->getVectorShape()(0),false);
892 : }
893 :
894 0 : return inputDummyDataVector_p;
895 : }
896 :
897 : // -----------------------------------------------------------------------
898 : //
899 : // -----------------------------------------------------------------------
900 0 : template<class T> Vector<T>& RegriddingKernel<T>::getOutputDataVector(DataCubeMap *outputData)
901 : {
902 0 : if (outputData->present(MS::DATA))
903 : {
904 0 : return outputData->getVector<T>(MS::DATA);
905 : }
906 0 : else if (not outputDummyDataVectorInitialized_p)
907 : {
908 0 : outputDummyDataVectorInitialized_p = true;
909 0 : outputDummyDataVector_p.resize(outputData->getVectorShape()(0),false);
910 : }
911 :
912 0 : return outputDummyDataVector_p;
913 : }
914 :
915 : //////////////////////////////////////////////////////////////////////////
916 : // DataInterpolationKernel class
917 : //////////////////////////////////////////////////////////////////////////
918 :
919 : // -----------------------------------------------------------------------
920 : //
921 : // -----------------------------------------------------------------------
922 0 : template<class T> DataInterpolationKernel<T>::DataInterpolationKernel( uInt interpolationMethod,
923 : Vector<Double> *inputFreq,
924 0 : Vector<Double> *outputFreq)
925 : {
926 0 : interpolationMethod_p = interpolationMethod;
927 0 : inputFreq_p = inputFreq;
928 0 : outputFreq_p = outputFreq;
929 0 : }
930 :
931 : // -----------------------------------------------------------------------
932 : //
933 : // -----------------------------------------------------------------------
934 0 : template<class T> void DataInterpolationKernel<T>::kernel( DataCubeMap *inputData,
935 : DataCubeMap *outputData)
936 : {
937 0 : Vector<T> &inputDataVector = getInputDataVector(inputData);
938 0 : Vector<T> &outputDataVector = getOutputDataVector(outputData);
939 :
940 0 : if (inputDataVector.size() > 1)
941 : {
942 0 : Vector<Bool> &inputFlagVector = getInputFlagVector(inputData);
943 0 : Vector<Bool> &outputFlagVector = getOutputFlagVector(outputData);
944 :
945 0 : InterpolateArray1D<Double,T>::interpolate( outputDataVector, // Output data
946 : outputFlagVector, // Output flags
947 0 : *outputFreq_p, // Out chan freq
948 0 : *inputFreq_p, // In chan freq
949 : inputDataVector, // Input data
950 : inputFlagVector, // Input Flags
951 0 : interpolationMethod_p, // Interpolation method
952 : false, // A good data point has its flag set to false
953 : false // If false extrapolated data points are set flagged
954 : );
955 : }
956 : else
957 : {
958 0 : outputDataVector = inputDataVector(0);
959 : }
960 :
961 0 : return;
962 : }
963 :
964 :
965 :
966 : //////////////////////////////////////////////////////////////////////////
967 : // DataFFTKernel class
968 : //////////////////////////////////////////////////////////////////////////
969 :
970 : // -----------------------------------------------------------------------
971 : //
972 : // -----------------------------------------------------------------------
973 0 : template<class T> DataFFTKernel<T>::DataFFTKernel(Double fftShift)
974 : {
975 0 : fftShift_p = fftShift;
976 0 : }
977 :
978 : // -----------------------------------------------------------------------
979 : //
980 : // -----------------------------------------------------------------------
981 0 : template<class T> void DataFFTKernel<T>::kernel(DataCubeMap *inputData,DataCubeMap *outputData)
982 : {
983 0 : Vector<T> &inputDataVector = getInputDataVector(inputData);
984 0 : Vector<T> &outputDataVector = getOutputDataVector(outputData);
985 :
986 0 : if (inputDataVector.size() > 1)
987 : {
988 0 : Vector<Bool> &inputFlagVector = getInputFlagVector(inputData);
989 0 : Vector<Bool> &outputFlagVector = getOutputFlagVector(outputData);
990 :
991 0 : fftshift(inputDataVector,inputFlagVector,outputDataVector,outputFlagVector);
992 : }
993 : else
994 : {
995 0 : outputDataVector = inputDataVector(0);
996 : }
997 :
998 0 : return;
999 : }
1000 :
1001 : // -----------------------------------------------------------------------
1002 : //
1003 : // -----------------------------------------------------------------------
1004 0 : template<class T> void DataFFTKernel<T>::fftshift( Vector<Complex> &inputDataVector,
1005 : Vector<Bool> &inputFlagVector,
1006 : Vector<Complex> &outputDataVector,
1007 : Vector<Bool> &outputFlagVector)
1008 : {
1009 0 : fFFTServer_p.fftshift( outputDataVector,
1010 : outputFlagVector,
1011 : (const Vector<T>)inputDataVector,
1012 : (const Vector<Bool>)inputFlagVector,
1013 : (const uInt)0, // In vectors axis 0 is the only dimension
1014 0 : (const Double)fftShift_p,
1015 : false, // A good data point has its flag set to false
1016 : false);
1017 :
1018 0 : return;
1019 : }
1020 :
1021 : // -----------------------------------------------------------------------
1022 : //
1023 : // -----------------------------------------------------------------------
1024 0 : template<class T> void DataFFTKernel<T>::fftshift( Vector<Float> &inputDataVector,
1025 : Vector<Bool> &inputFlagVector,
1026 : Vector<Float> &outputDataVector,
1027 : Vector<Bool> &outputFlagVector)
1028 : {
1029 0 : fFFTServer_p.fftshift( outputDataVector,
1030 : outputFlagVector,
1031 : (const Vector<T>)inputDataVector,
1032 : (const Vector<Bool>)inputFlagVector,
1033 : (const uInt)0, // In vectors axis 0 is the only dimension
1034 0 : (const Double)fftShift_p,
1035 : false); // A good data point has its flag set to false
1036 :
1037 0 : return;
1038 : }
1039 :
1040 : } //# NAMESPACE VI - END
1041 :
1042 : } //# NAMESPACE CASA - END
1043 :
1044 :
|
