Line data Source code
1 : //# ImagePolarimetry.cc: polarimetric analysis
2 : //# Copyright (C) 1996,1997,1998,1999,2000,2001,2002,2003
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: ImagePolarimetry.cc 20652 2009-07-06 05:04:32Z Malte.Marquarding $
27 :
28 : #include <casacore/casa/OS/Timer.h>
29 :
30 : #include <imageanalysis/ImageAnalysis/ImagePolarimetry.h>
31 :
32 : #include <casacore/casa/Arrays/Array.h>
33 : #include <casacore/casa/Arrays/ArrayMath.h>
34 : #include <casacore/casa/Arrays/Vector.h>
35 : #include <casacore/casa/Arrays/Matrix.h>
36 : #include <casacore/casa/Arrays/MaskedArray.h>
37 : #include <casacore/casa/Arrays/MaskArrMath.h>
38 : #include <casacore/coordinates/Coordinates/CoordinateSystem.h>
39 : #include <casacore/coordinates/Coordinates/StokesCoordinate.h>
40 : #include <casacore/coordinates/Coordinates/LinearCoordinate.h>
41 : #include <casacore/casa/Exceptions/Error.h>
42 : #include <casacore/scimath/Functionals/Polynomial.h>
43 : #include <casacore/images/Images/ImageInterface.h>
44 : #include <casacore/images/Images/SubImage.h>
45 : #include <casacore/images/Images/ImageExpr.h>
46 : #include <imageanalysis/ImageAnalysis/ImageFFT.h>
47 : #include <casacore/images/Regions/ImageRegion.h>
48 : #include <casacore/images/Images/ImageSummary.h>
49 : #include <casacore/images/Images/TempImage.h>
50 : #include <casacore/lattices/Lattices/Lattice.h>
51 : #include <casacore/lattices/LRegions/LCSlicer.h>
52 : #include <casacore/lattices/LEL/LatticeExprNode.h>
53 : #include <casacore/lattices/LEL/LatticeExpr.h>
54 : #include <casacore/lattices/Lattices/TiledLineStepper.h>
55 : #include <casacore/lattices/Lattices/LatticeStepper.h>
56 : #include <casacore/lattices/Lattices/LatticeIterator.h>
57 : #include <casacore/lattices/Lattices/LatticeUtilities.h>
58 : #include <casacore/lattices/Lattices/MaskedLatticeIterator.h>
59 : #include <casacore/lattices/LatticeMath/LatticeStatistics.h>
60 : #include <casacore/lattices/LRegions/LCPagedMask.h>
61 : #include <casacore/casa/Logging/LogIO.h>
62 : #include <casacore/casa/Logging/LogOrigin.h>
63 : #include <casacore/casa/BasicMath/Math.h>
64 : #include <casacore/casa/BasicSL/Constants.h>
65 : #include <casacore/scimath/Mathematics/NumericTraits.h>
66 : #include <casacore/casa/System/ProgressMeter.h>
67 : #include <casacore/casa/Quanta/QC.h>
68 : #include <casacore/casa/Quanta/MVAngle.h>
69 : #include <casacore/casa/Utilities/GenSort.h>
70 : #include <casacore/casa/Utilities/Assert.h>
71 : #include <casacore/casa/BasicSL/String.h>
72 :
73 : #include <sstream>
74 :
75 : using namespace casacore;
76 : namespace casa {
77 :
78 : const std::map<ImagePolarimetry::StokesTypes, String> ImagePolarimetry::polMap = {
79 : {I, "I"}, {Q, "Q"}, {U, "U"}, {V, "V"}
80 : };
81 :
82 44 : ImagePolarimetry::ImagePolarimetry (const ImageInterface<Float>& image)
83 68 : : _image(image.cloneII())
84 : {
85 44 : _stokes.resize(4);
86 44 : _stokesStats.resize(4);
87 44 : _stokes.set(0);
88 44 : _stokesStats.set(0);
89 44 : _findStokes();
90 38 : _createBeamsEqMat();
91 38 : }
92 :
93 0 : ImagePolarimetry::ImagePolarimetry(const ImagePolarimetry &other) {
94 0 : operator=(other);
95 0 : }
96 :
97 0 : ImagePolarimetry &ImagePolarimetry::operator=(const ImagePolarimetry &other) {
98 0 : if (this != &other) {
99 0 : _image.reset(other._image->cloneII());
100 0 : const auto n = _stokes.size();
101 0 : for (size_t i=0; i<n; ++i) {
102 0 : if (_stokes[i]) {
103 0 : delete _stokes[i];
104 0 : _stokes[i] = nullptr;
105 : }
106 0 : if (other._stokes[i]) {
107 0 : _stokes[i] = other._stokes[i]->cloneII();
108 : }
109 : }
110 : // Just delete fitter. It will make a new one when needed.
111 0 : if (_fitter) {
112 0 : delete _fitter;
113 0 : _fitter = nullptr;
114 : }
115 : // Remake Statistics objects as needed
116 0 : _oldClip = 0.0;
117 0 : for (size_t i=0; i<n; ++i) {
118 0 : if (_stokesStats[i]) {
119 0 : delete _stokesStats[i];
120 0 : _stokesStats[i] = nullptr;
121 : }
122 : }
123 0 : _beamsEqMat.assign(other._beamsEqMat);
124 : }
125 0 : return *this;
126 : }
127 :
128 38 : ImagePolarimetry::~ImagePolarimetry() {
129 38 : _cleanup();
130 38 : }
131 :
132 0 : ImageExpr<Complex> ImagePolarimetry::complexLinearPolarization() {
133 0 : _hasQU();
134 0 : _checkQUBeams(false);
135 : LatticeExprNode node(
136 : casacore::formComplex(
137 0 : *_stokes[ImagePolarimetry::Q], *_stokes[ImagePolarimetry::U]
138 : )
139 0 : );
140 0 : LatticeExpr<Complex> le(node);
141 0 : ImageExpr<Complex> ie(le, String("ComplexLinearPolarization"));
142 : // Need a Complex Linear Polarization type
143 0 : _fiddleStokesCoordinate(ie, Stokes::Plinear);
144 0 : ie.setUnits(_image->units());
145 0 : _setInfo(ie, Q);
146 0 : return ie;
147 : }
148 :
149 0 : void ImagePolarimetry::_setInfo(
150 : ImageInterface<Complex>& im, StokesTypes stokes
151 : ) const {
152 0 : auto info = _image->imageInfo();
153 0 : if (info.hasMultipleBeams()) {
154 0 : info.setBeams(_stokes[stokes]->imageInfo().getBeamSet());
155 : }
156 0 : im.setImageInfo(info);
157 0 : }
158 :
159 0 : void ImagePolarimetry::_setInfo(
160 : ImageInterface<Float>& im, StokesTypes stokes
161 : ) const {
162 0 : auto info = _image->imageInfo();
163 0 : if (info.hasMultipleBeams()) {
164 0 : info.setBeams(_stokes[stokes]->imageInfo().getBeamSet());
165 : }
166 0 : im.setImageInfo(info);
167 0 : }
168 :
169 0 : ImageExpr<Complex> ImagePolarimetry::complexFractionalLinearPolarization() {
170 0 : LogIO os(LogOrigin("ImagePolarimetry", __func__, WHERE));
171 0 : _hasQU();
172 0 : ThrowIf(
173 : ! _stokes[ImagePolarimetry::I],
174 : "This image does not have Stokes I so cannot "
175 : "provide fractional linear polarization"
176 : );
177 0 : _checkIQUBeams(false);
178 : LatticeExprNode nodeQU(
179 : casacore::formComplex(
180 0 : *_stokes[ImagePolarimetry::Q], *_stokes[ImagePolarimetry::U]
181 : )
182 0 : );
183 0 : LatticeExprNode nodeI(*_stokes[ImagePolarimetry::I]);
184 0 : LatticeExpr<Complex> le(nodeQU/nodeI);
185 0 : ImageExpr<Complex> ie(le, String("ComplexFractionalLinearPolarization"));
186 : // Need a Complex Linear Polarization type
187 0 : _fiddleStokesCoordinate(ie, Stokes::PFlinear);
188 0 : ie.setUnits(Unit(""));
189 0 : _setInfo(ie, I);
190 0 : return ie;
191 : }
192 :
193 0 : ImageExpr<Float> ImagePolarimetry::fracLinPol(
194 : Bool debias, Float clip, Float sigma
195 : ) {
196 0 : LogIO os(LogOrigin("ImagePolarimetry", __func__, WHERE));
197 0 : _hasQU();
198 0 : ThrowIf(
199 : ! _stokes[ImagePolarimetry::I],
200 : "This image does not have Stokes I so cannot "
201 : "provide fractional linear polarization"
202 : );
203 0 : Vector<StokesTypes> types(3);
204 0 : types[0] = I; types[1] = Q; types[2] = U;
205 0 : _checkIQUBeams(false);
206 0 : auto nodePol = _makePolIntNode(os, debias, clip, sigma, true, false);
207 0 : LatticeExprNode nodeI(*_stokes[ImagePolarimetry::I]);
208 0 : LatticeExpr<Float> le(nodePol/nodeI);
209 0 : ImageExpr<Float> ie(le, String("FractionalLinearPolarization"));
210 0 : ie.setUnits(Unit(""));
211 0 : auto ii = _image->imageInfo();
212 0 : ii.removeRestoringBeam();
213 0 : ie.setImageInfo(ii);
214 0 : _fiddleStokesCoordinate(ie, Stokes::PFlinear);
215 0 : return ie;
216 : }
217 :
218 0 : ImageExpr<Float> ImagePolarimetry::sigmaFracLinPol(Float clip, Float sigma) {
219 : // sigma_m = m * sqrt( (sigmaP/p)**2 + (sigmaI/I)**2) )
220 : // sigmaP = sigmaQU
221 : // sigmaI = sigmaI
222 0 : LogIO os(LogOrigin("ImagePolarimetry", __func__, WHERE));
223 0 : _hasQU();
224 0 : ThrowIf(
225 : ! _stokes[ImagePolarimetry::I],
226 : "This image does not have Stokes I so cannot provide "
227 : "fractional linear polarization"
228 : );
229 0 : _checkIQUBeams(false);
230 : // Make nodes. Don't bother debiasing.
231 0 : Bool debias = false;
232 0 : auto nodePol = _makePolIntNode(os, debias, clip, sigma, true, false);
233 0 : LatticeExprNode nodeI(*_stokes[ImagePolarimetry::I]);
234 :
235 : // Make expression. We assume sigmaI = sigmaQU which is true with
236 : // no dynamic range limititation. Perhaps we should work out
237 : // sigmaI as well.
238 :
239 0 : const auto sigma2 = sigmaLinPolInt(clip, sigma);
240 0 : LatticeExprNode n0(nodePol/nodeI);
241 0 : LatticeExprNode n1(pow(sigma2/nodePol,2));
242 0 : LatticeExprNode n2(pow(sigma2/nodeI,2));
243 0 : LatticeExpr<Float> le(n0 * sqrt(n1 + n2));
244 0 : ImageExpr<Float> ie(le, String("FractionalLinearPolarizationError"));
245 0 : ie.setUnits(Unit(""));
246 0 : auto ii = _image->imageInfo();
247 0 : ii.removeRestoringBeam();
248 0 : ie.setImageInfo(ii);
249 0 : _fiddleStokesCoordinate(ie, Stokes::PFlinear);
250 0 : return ie;
251 : }
252 :
253 0 : ImageExpr<Float> ImagePolarimetry::fracTotPol(
254 : Bool debias, Float clip, Float sigma
255 : ) {
256 0 : LogIO os(LogOrigin("ImagePolarimetry", __func__, WHERE));
257 : Bool doLin, doCirc;
258 0 : _setDoLinDoCirc(doLin, doCirc);
259 0 : auto nodePol = _makePolIntNode(os, debias, clip, sigma, doLin, doCirc);
260 0 : LatticeExprNode nodeI(*_stokes[ImagePolarimetry::I]);
261 0 : LatticeExpr<Float> le(nodePol/nodeI);
262 0 : ImageExpr<Float> ie(le, String("FractionalTotalPolarization"));
263 0 : ie.setUnits(Unit(""));
264 0 : auto ii = _image->imageInfo();
265 0 : ii.removeRestoringBeam();
266 0 : ie.setImageInfo(ii);
267 0 : _fiddleStokesCoordinate(ie, Stokes::PFtotal);
268 0 : return ie;
269 : }
270 :
271 0 : void ImagePolarimetry::_setDoLinDoCirc(Bool& doLin, Bool& doCirc) const {
272 0 : LogIO os(LogOrigin("ImagePolarimetry", __func__, WHERE));
273 0 : doLin = _stokes[ImagePolarimetry::Q] && _stokes[ImagePolarimetry::U];
274 0 : doCirc = _stokes[ImagePolarimetry::V];
275 : // Should never happen
276 0 : AlwaysAssert((doLin || doCirc), AipsError);
277 0 : ThrowIf(
278 : ! _stokes[ImagePolarimetry::I],
279 : "This image does not have Stokes I so this calculation "
280 : "cannot be carried out"
281 : );
282 0 : if (doLin && ! _checkIQUBeams(false, false)) {
283 : os << LogIO::WARN << "I, Q, and U beams are not the same, cannot do "
284 0 : << "linear portion" << LogIO::POST;
285 0 : doLin = false;
286 : }
287 0 : if (doCirc && ! _checkIVBeams(false, false)) {
288 : os << LogIO::WARN << "I and V beams are not the same, cannot do "
289 0 : << "circular portion" << LogIO::POST;
290 0 : doCirc = false;
291 : }
292 0 : ThrowIf(
293 : ! (doLin || doCirc), "Can do neither linear nor circular portions"
294 : );
295 0 : }
296 :
297 0 : ImageExpr<Float> ImagePolarimetry::sigmaFracTotPol(Float clip, Float sigma) {
298 : // sigma_m = m * sqrt( (sigmaP/P)**2 + (sigmaI/I)**2) )
299 : // sigmaP = sigmaQU
300 : // sigmaI = sigmaI
301 0 : LogIO os(LogOrigin("ImagePolarimetry", __func__, WHERE));
302 : Bool doLin, doCirc;
303 0 : _setDoLinDoCirc(doLin, doCirc);
304 : // Make nodes. Don't bother debiasing.
305 0 : Bool debias = false;
306 0 : auto nodePol = _makePolIntNode(os, debias, clip, sigma, doLin, doCirc);
307 0 : LatticeExprNode nodeI(*_stokes[ImagePolarimetry::I]);
308 : // Make expression. We assume sigmaI = sigmaQU which is true with
309 : // no dynamic range limitation. Perhaps we should work out
310 : // sigmaI as well.
311 0 : const auto sigma2 = sigmaTotPolInt(clip, sigma);
312 0 : LatticeExprNode n0(nodePol/nodeI);
313 0 : LatticeExprNode n1(pow(sigma2/nodePol,2));
314 0 : LatticeExprNode n2(pow(sigma2/nodeI,2));
315 0 : LatticeExpr<Float> le(n0 * sqrt(n1 + n2));
316 0 : ImageExpr<Float> ie(le, String("FractionalLinearPolarizationError"));
317 0 : ie.setUnits(Unit(""));
318 0 : auto ii = _image->imageInfo();
319 0 : ii.removeRestoringBeam();
320 0 : ie.setImageInfo(ii);
321 0 : _fiddleStokesCoordinate(ie, Stokes::PFlinear);
322 0 : return ie;
323 : }
324 :
325 0 : void ImagePolarimetry::fourierRotationMeasure(
326 : ImageInterface<Complex>& cpol, Bool zeroZeroLag
327 : ) {
328 0 : LogIO os(LogOrigin("ImagePolarimetry", __func__, WHERE));
329 0 : _hasQU();
330 0 : _checkQUBeams(true, false);
331 0 : CoordinateSystem dCS;
332 0 : Stokes::StokesTypes dType = Stokes::Plinear;
333 0 : const auto shape = singleStokesShape(dCS, dType);
334 0 : if(! cpol.shape().isEqual(shape)) {
335 0 : os << "The provided image has the wrong shape " << cpol.shape()
336 0 : << endl;
337 0 : os << "It should be of shape " << shape << LogIO::EXCEPTION;
338 : }
339 0 : const auto& cSys = _image->coordinates();
340 : Int spectralCoord, fAxis;
341 0 : _findFrequencyAxis (spectralCoord, fAxis, cSys, -1);
342 : // Make Complex (Q,U) image
343 0 : LatticeExprNode node;
344 0 : if (zeroZeroLag) {
345 : TempImage<Float> tQ(
346 0 : _stokes[ImagePolarimetry::Q]->shape(),
347 0 : _stokes[ImagePolarimetry::Q]->coordinates()
348 0 : );
349 0 : if (_stokes[ImagePolarimetry::Q]->isMasked()) {
350 0 : tQ.makeMask(String("mask0"), true, true, false, false);
351 : }
352 0 : LatticeUtilities::copyDataAndMask(
353 0 : os, tQ, *_stokes[ImagePolarimetry::Q], false
354 : );
355 0 : _subtractProfileMean (tQ, fAxis);
356 : TempImage<Float> tU(
357 0 : _stokes[ImagePolarimetry::U]->shape(),
358 0 : _stokes[ImagePolarimetry::U]->coordinates()
359 0 : );
360 0 : if (_stokes[ImagePolarimetry::U]->isMasked()) {
361 0 : tU.makeMask(String("mask0"), true, true, false, false);
362 : }
363 0 : LatticeUtilities::copyDataAndMask(
364 0 : os, tU, *_stokes[ImagePolarimetry::U], false
365 : );
366 0 : _subtractProfileMean (tU, fAxis);
367 : // The TempImages will be cloned be LatticeExprNode so it's ok
368 : // that they go out of scope
369 0 : node = LatticeExprNode(formComplex(tQ, tU));
370 : }
371 : else {
372 0 : node = LatticeExprNode(
373 : formComplex(
374 0 : *_stokes[ImagePolarimetry::Q], *_stokes[ImagePolarimetry::U]
375 : )
376 0 : );
377 : }
378 0 : LatticeExpr<Complex> le(node);
379 0 : ImageExpr<Complex> ie(le, String("ComplexLinearPolarization"));
380 : // Do FFT of spectral coordinate
381 0 : Vector<Bool> axes(ie.ndim(),false);
382 0 : axes(fAxis) = true;
383 0 : ImageFFT<Complex> fftserver;
384 0 : fftserver.fft(ie, axes);
385 : // Recover Complex result. Coordinates are updated to include Fourier
386 : // coordinate, miscellaneous things (MiscInfo, ImageInfo, units, history)
387 : // and mask (if output has one) are copied to cpol
388 0 : fftserver.getComplex(cpol);
389 : // Fiddle time coordinate to be a RotationMeasure coordinate
390 : auto f = _findCentralFrequency(
391 0 : cSys.coordinate(spectralCoord), ie.shape()(fAxis)
392 0 : );
393 0 : _fiddleTimeCoordinate(cpol, f, spectralCoord);
394 : // Set Stokes coordinate to be correct type
395 0 : _fiddleStokesCoordinate(cpol, Stokes::Plinear);
396 : // Set units and ImageInfo
397 0 : cpol.setUnits(_image->units());
398 0 : _setInfo(cpol, Q);
399 0 : }
400 :
401 0 : Float ImagePolarimetry::sigmaLinPolInt(Float clip, Float sigma) {
402 : // sigma_P = sigma_QU
403 0 : LogIO os(LogOrigin("ImagePolarimetry", __func__, WHERE));
404 0 : ThrowIf(
405 : ! _stokes[ImagePolarimetry::Q] && ! _stokes[ImagePolarimetry::U],
406 : "This image does not have Stokes Q and U so cannot "
407 : "provide linear polarization"
408 : );
409 0 : _checkQUBeams(false);
410 0 : Float sigma2 = 0.0;
411 0 : if (sigma > 0) {
412 0 : sigma2 = sigma;
413 : }
414 : else {
415 0 : os << LogIO::NORMAL << "Determined noise from Q&U images to be ";
416 0 : auto sq = _sigma(ImagePolarimetry::Q, clip);
417 0 : auto su = _sigma(ImagePolarimetry::U, clip);
418 0 : sigma2 = (sq+su)/2.0;
419 : }
420 0 : return sigma2;
421 : }
422 :
423 16 : ImageExpr<Float> ImagePolarimetry::linPolPosAng(Bool radians) const {
424 48 : LogIO os(LogOrigin("ImagePolarimetry", __func__, WHERE));
425 16 : ThrowIf(
426 : ! _stokes[ImagePolarimetry::Q] && ! _stokes[ImagePolarimetry::U],
427 : "This image does not have Stokes Q and U so cannot "
428 : "provide linear polarization"
429 : );
430 16 : _checkQUBeams(false);
431 : // Make expression. LEL function "pa" returns degrees
432 16 : Float fac = radians ? C::pi / 180.0 : 1.0;
433 : LatticeExprNode node(
434 32 : fac*pa(*_stokes[ImagePolarimetry::U], *_stokes[ImagePolarimetry::Q])
435 48 : );
436 32 : LatticeExpr<Float> le(node);
437 32 : ImageExpr<Float> ie(le, String("LinearlyPolarizedPositionAngle"));
438 16 : ie.setUnits(Unit(radians ? "rad" : "deg"));
439 32 : auto ii = _image->imageInfo();
440 16 : ii.removeRestoringBeam();
441 16 : ie.setImageInfo(ii);
442 16 : _fiddleStokesCoordinate(ie, Stokes::Pangle);
443 32 : return ie;
444 : }
445 :
446 16 : ImageExpr<Float> ImagePolarimetry::sigmaLinPolPosAng(
447 : Bool radians, Float clip, Float sigma
448 : ) {
449 : // sigma_PA = sigmaQU / 2P
450 16 : ThrowIf(
451 : ! (_stokes[ImagePolarimetry::Q] || _stokes[ImagePolarimetry::U]),
452 : "This image does not have Stokes Q and U so "
453 : "cannot provide linear polarization"
454 : );
455 16 : _checkQUBeams(false);
456 16 : Float sigma2 = sigma > 0 ? sigma : this->sigma(clip);
457 16 : Float fac = 0.5 * sigma2;
458 16 : if (! radians) {
459 0 : fac *= 180 / C::pi;
460 : }
461 : LatticeExprNode node(
462 32 : fac / amp(*_stokes[ImagePolarimetry::U], *_stokes[ImagePolarimetry::Q])
463 48 : );
464 32 : LatticeExpr<Float> le(node);
465 32 : ImageExpr<Float> ie(le, String("LinearlyPolarizedPositionAngleError"));
466 16 : ie.setUnits(Unit(radians ? "rad" : "deg"));
467 32 : auto ii = _image->imageInfo();
468 16 : ii.removeRestoringBeam();
469 16 : ie.setImageInfo(ii);
470 16 : _fiddleStokesCoordinate(ie, Stokes::Pangle);
471 32 : return ie;
472 : }
473 :
474 12 : Float ImagePolarimetry::sigma(Float clip) {
475 24 : LogIO os(LogOrigin("ImagePolarimetry", __func__, WHERE));
476 12 : Float sigma2 = 0.0;
477 12 : if (_stokes[ImagePolarimetry::V]) {
478 12 : os << LogIO::NORMAL << "Determined noise from V image to be ";
479 12 : sigma2 = _sigma(ImagePolarimetry::V, clip);
480 : }
481 0 : else if (
482 0 : _stokes[ImagePolarimetry::Q] && _stokes[ImagePolarimetry::U]
483 0 : && _checkQUBeams(false, false)
484 : ) {
485 0 : sigma2 = sigmaLinPolInt(clip);
486 : }
487 0 : else if (_stokes[ImagePolarimetry::Q]) {
488 : os << LogIO::NORMAL << "Determined noise from Q image to be "
489 0 : << LogIO::POST;
490 0 : sigma2 = _sigma(ImagePolarimetry::Q, clip);
491 : }
492 0 : else if (_stokes[ImagePolarimetry::U]) {
493 : os << LogIO::NORMAL << "Determined noise from U image to be "
494 0 : << LogIO::POST;
495 0 : sigma2 = _sigma(ImagePolarimetry::U, clip);
496 : }
497 0 : else if (_stokes[ImagePolarimetry::I]!=0) {
498 : os << LogIO::NORMAL << "Determined noise from I image to be "
499 0 : << LogIO::POST;
500 0 : sigma2 = _sigma(ImagePolarimetry::I, clip);
501 : }
502 12 : os << sigma2 << LogIO::POST;
503 24 : return sigma2;
504 : }
505 :
506 16 : void ImagePolarimetry::rotationMeasure(
507 : ImageInterface<Float>*& rmOutPtr, ImageInterface<Float>*& rmOutErrorPtr,
508 : ImageInterface<Float>*& pa0OutPtr, ImageInterface<Float>*& pa0OutErrorPtr,
509 : ImageInterface<Float>*& nTurnsOutPtr, ImageInterface<Float>*& chiSqOutPtr,
510 : Int axis, Float rmMax, Float maxPaErr, Float sigma, Float rmFg,
511 : Bool showProgress
512 : ) {
513 48 : LogIO os(LogOrigin("ImagePolarimetry", __func__, WHERE));
514 16 : _hasQU();
515 16 : _checkQUBeams(false);
516 : // Do we have anything to do ?
517 16 : ThrowIf(
518 : ! (rmOutPtr || rmOutErrorPtr || pa0OutPtr || pa0OutErrorPtr),
519 : "No output images specified"
520 : );
521 : // Find expected shape of output RM images (Stokes and spectral axes gone)
522 32 : CoordinateSystem cSysRM;
523 : Int fAxis, sAxis;
524 32 : const auto shapeRM = rotationMeasureShape(cSysRM, fAxis, sAxis, os, axis);
525 32 : const auto shapeNTurns = shapeRM;
526 32 : const auto shapeChiSq = shapeRM;
527 : // Check RM image shapes
528 16 : if (rmOutPtr && ! rmOutPtr->shape().isEqual(shapeRM)) {
529 : os << "The provided Rotation Measure image has the wrong shape "
530 0 : << rmOutPtr->shape() << endl;
531 0 : os << "It should be of shape " << shapeRM << LogIO::EXCEPTION;
532 : }
533 16 : if (rmOutErrorPtr && !rmOutErrorPtr->shape().isEqual(shapeRM)) {
534 : os << "The provided Rotation Measure error image has the wrong shape "
535 0 : << rmOutErrorPtr->shape() << endl;
536 0 : os << "It should be of shape " << shapeRM << LogIO::EXCEPTION;
537 : }
538 : // Check position angle image shapes
539 32 : CoordinateSystem cSysPA;
540 32 : const auto shapePA = positionAngleShape(cSysPA, fAxis, sAxis, os, axis);
541 16 : if (pa0OutPtr && ! pa0OutPtr->shape().isEqual(shapePA)) {
542 : os << "The provided position angle at zero wavelength image has the "
543 0 : << "wrong shape " << pa0OutPtr->shape() << endl;
544 0 : os << "It should be of shape " << shapePA << LogIO::EXCEPTION;
545 : }
546 16 : if (pa0OutErrorPtr && ! pa0OutErrorPtr->shape().isEqual(shapePA)) {
547 : os << "The provided position angle at zero wavelength image has the "
548 0 : << "wrong shape " << pa0OutErrorPtr->shape() << endl;
549 0 : os << "It should be of shape " << shapePA << LogIO::EXCEPTION;
550 : }
551 16 : if (nTurnsOutPtr && ! nTurnsOutPtr->shape().isEqual(shapeNTurns)) {
552 : os << "The provided nTurns image has the wrong shape "
553 0 : << nTurnsOutPtr->shape() << endl;
554 0 : os << "It should be of shape " << shapeNTurns << LogIO::EXCEPTION;
555 : }
556 16 : if (chiSqOutPtr && ! chiSqOutPtr->shape().isEqual(shapeChiSq)) {
557 : os << "The provided chi squared image has the wrong shape "
558 0 : << chiSqOutPtr->shape() << endl;
559 0 : os << "It should be of shape " << shapeChiSq << LogIO::EXCEPTION;
560 : }
561 : // Generate linear polarization position angle image expressions
562 : // and error in radians
563 16 : Bool radians = true;
564 16 : Float clip = 10.0;
565 32 : const auto pa = linPolPosAng(radians);
566 32 : const auto paerr = sigmaLinPolPosAng(radians, clip, sigma);
567 32 : CoordinateSystem cSys0 = pa.coordinates();
568 : // Set frequency axis units to Hz
569 16 : auto fAxisWorld = cSys0.pixelAxisToWorldAxis(fAxis);
570 16 : ThrowIf(
571 : fAxisWorld < 0,
572 : "World axis has been removed for the frequency pixel axis"
573 : );
574 32 : auto axisUnits = cSys0.worldAxisUnits();
575 16 : axisUnits(fAxisWorld) = String("Hz");
576 16 : ThrowIf(
577 : ! cSys0.setWorldAxisUnits(axisUnits),
578 : "Failed to set frequency axis units to Hz because "
579 : + cSys0.errorMessage()
580 : );
581 : // Do we have enough frequency pixels ?
582 16 : const uInt nFreq = pa.shape()(fAxis);
583 19 : ThrowIf(
584 : nFreq < 3,
585 : "This image only has " + String::toString(nFreq)
586 : + " frequencies, this is not enough"
587 : );
588 : // Copy units only over. The output images don't have a beam
589 : // so unset beam. MiscInfo and history require writable II.
590 : // We leave this to the caller who knows what sort of II these are.
591 30 : auto ii = _image->imageInfo();
592 15 : ii.removeRestoringBeam();
593 15 : if (rmOutPtr) {
594 15 : rmOutPtr->setImageInfo(ii);
595 15 : rmOutPtr->setUnits(Unit("rad/m/m"));
596 : }
597 15 : if (rmOutErrorPtr) {
598 0 : rmOutErrorPtr->setImageInfo(ii);
599 0 : rmOutErrorPtr->setUnits(Unit("rad/m/m"));
600 : }
601 15 : if (pa0OutPtr) {
602 1 : pa0OutPtr->setImageInfo(ii);
603 1 : pa0OutPtr->setUnits(Unit("deg"));
604 : }
605 15 : if (pa0OutErrorPtr) {
606 0 : pa0OutErrorPtr->setImageInfo(ii);
607 0 : pa0OutErrorPtr->setUnits(Unit("deg"));
608 : }
609 15 : if (nTurnsOutPtr) {
610 0 : nTurnsOutPtr->setImageInfo(ii);
611 0 : nTurnsOutPtr->setUnits(Unit(""));
612 : }
613 15 : if (chiSqOutPtr) {
614 0 : chiSqOutPtr->setImageInfo(ii);
615 0 : chiSqOutPtr->setUnits(Unit(""));
616 : }
617 : // Get lambda squared in m**2
618 30 : Vector<Double> freqs(nFreq);
619 30 : Vector<Float> wsq(nFreq);
620 30 : Vector<Double> world;
621 45 : Vector<Double> pixel(cSys0.referencePixel().copy());
622 15 : Double c = QC::c( ).getValue(Unit("m/s"));
623 15 : Double csq = c*c;
624 360 : for (uInt i=0; i<nFreq; ++i) {
625 345 : pixel(fAxis) = i;
626 345 : ThrowIf(
627 : !cSys0.toWorld(world, pixel),
628 : "Failed to convert pixel to world because "
629 : + cSys0.errorMessage()
630 : );
631 345 : freqs(i) = world(fAxisWorld);
632 : // m**2
633 345 : wsq(i) = csq / freqs(i) / freqs(i);
634 : }
635 : // Sort into increasing wavelength
636 30 : Vector<uInt> sortidx;
637 15 : GenSortIndirect<Float>::sort(
638 : sortidx, wsq, Sort::Ascending, Sort::QuickSort|Sort::NoDuplicates
639 : );
640 30 : Vector<Float> wsqsort(sortidx.size());
641 360 : for (uInt i=0; i<wsqsort.size(); ++i) {
642 345 : wsqsort[i] = wsq[sortidx[i]];
643 : }
644 : // Make fitter
645 15 : if (! _fitter) {
646 15 : _fitter = new LinearFitSVD<Float>;
647 : // Create and set the polynomial functional
648 : // p = c(0) + c(1)*x where x = lambda**2
649 : // PA = PA0 + RM*Lambda**2
650 30 : Polynomial<AutoDiff<Float> > poly1(1);
651 : // Makes a copy of poly1
652 15 : _fitter->setFunction(poly1);
653 : }
654 : // Deal with masks. The outputs are all given a mask if possible as we
655 : // don't know at this point whether output points will be masked or not
656 30 : IPosition whereRM;
657 15 : auto isMaskedRM = false;
658 15 : Lattice<Bool>* outRMMaskPtr = nullptr;
659 15 : if (rmOutPtr) {
660 15 : isMaskedRM = _dealWithMask (outRMMaskPtr, rmOutPtr, os, "RM");
661 15 : whereRM.resize(rmOutPtr->ndim());
662 15 : whereRM = 0;
663 : }
664 15 : auto isMaskedRMErr = false;
665 15 : Lattice<Bool>* outRMErrMaskPtr = nullptr;
666 15 : if (rmOutErrorPtr) {
667 0 : isMaskedRMErr = _dealWithMask(
668 0 : outRMErrMaskPtr, rmOutErrorPtr, os, String("RM error")
669 : );
670 0 : whereRM.resize(rmOutErrorPtr->ndim());
671 0 : whereRM = 0;
672 : }
673 30 : IPosition wherePA;
674 15 : auto isMaskedPa0 = false;
675 15 : Lattice<Bool>* outPa0MaskPtr = nullptr;
676 15 : if (pa0OutPtr) {
677 1 : isMaskedPa0 = _dealWithMask(
678 2 : outPa0MaskPtr, pa0OutPtr, os, String("Position Angle")
679 : );
680 1 : wherePA.resize(pa0OutPtr->ndim());
681 1 : wherePA = 0;
682 : }
683 15 : auto isMaskedPa0Err = false;
684 15 : Lattice<Bool>* outPa0ErrMaskPtr = nullptr;
685 15 : if (pa0OutErrorPtr) {
686 0 : isMaskedPa0Err = _dealWithMask(
687 : outPa0ErrMaskPtr, pa0OutErrorPtr, os, "Position Angle error"
688 : );
689 0 : wherePA.resize(pa0OutErrorPtr->ndim());
690 0 : wherePA = 0;
691 : }
692 30 : IPosition whereNTurns;
693 15 : auto isMaskedNTurns = false;
694 15 : Lattice<Bool>* outNTurnsMaskPtr = 0;
695 15 : if (nTurnsOutPtr) {
696 0 : isMaskedNTurns = _dealWithMask(
697 : outNTurnsMaskPtr, nTurnsOutPtr, os, "nTurns"
698 : );
699 0 : whereNTurns.resize(nTurnsOutPtr->ndim());
700 0 : whereNTurns = 0;
701 : }
702 30 : IPosition whereChiSq;
703 15 : auto isMaskedChiSq = false;
704 15 : Lattice<Bool>* outChiSqMaskPtr = nullptr;
705 15 : if (chiSqOutPtr) {
706 0 : isMaskedChiSq = _dealWithMask(
707 0 : outChiSqMaskPtr, chiSqOutPtr, os, String("chi sqared")
708 : );
709 0 : whereChiSq.resize(chiSqOutPtr->ndim());
710 0 : whereChiSq = 0;
711 : }
712 30 : Array<Bool> tmpMaskRM(IPosition(shapeRM.size(), 1), true);
713 30 : Array<Float> tmpValueRM(IPosition(shapeRM.size(), 1), 0.0f);
714 30 : Array<Bool> tmpMaskPA(IPosition(shapePA.size(), 1), true);
715 30 : Array<Float> tmpValuePA(IPosition(shapePA.size(), 1), 0.0f);
716 30 : Array<Float> tmpValueNTurns(IPosition(shapeNTurns.size(), 1), 0.0f);
717 30 : Array<Bool> tmpMaskNTurns(IPosition(shapeNTurns.size(), 1), true);
718 30 : Array<Float> tmpValueChiSq(IPosition(shapeChiSq.size(), 1), 0.0f);
719 30 : Array<Bool> tmpMaskChiSq(IPosition(shapeChiSq.size(), 1), true);
720 : // Iterate
721 30 : const IPosition tileShape = pa.niceCursorShape();
722 30 : TiledLineStepper ts(pa.shape(), tileShape, fAxis);
723 30 : RO_MaskedLatticeIterator<Float> it(pa, ts);
724 : Float rm, rmErr, pa0, pa0Err, rChiSq, nTurns;
725 : uInt j, k, l, m;
726 15 : static const Double degPerRad = 180/C::pi;
727 15 : maxPaErr /= degPerRad;
728 15 : maxPaErr = abs(maxPaErr);
729 15 : Bool doRM = whereRM.size() > 0;
730 15 : Bool doPA = wherePA.size() > 0;
731 15 : Bool doNTurns = whereNTurns.size() > 0;
732 15 : Bool doChiSq = whereChiSq.size() > 0;
733 15 : unique_ptr<ProgressMeter> pProgressMeter;
734 15 : if (showProgress) {
735 15 : Double nMin = 0.0;
736 15 : Double nMax = 1.0;
737 75 : for (Int i=0; i<Int(pa.ndim()); ++i) {
738 60 : if (i!=fAxis) {
739 45 : nMax *= pa.shape()[i];
740 : }
741 : }
742 30 : pProgressMeter.reset(
743 : new ProgressMeter(
744 : nMin, nMax, "Profiles fitted", "Fitting",
745 15 : "", "", true, max(1,Int(nMax/100))
746 15 : )
747 : );
748 : }
749 : // As a (temporary?) workaround the cache of the main image is set up in
750 : // such a way that it can hold the full frequency and stokes axes.
751 : // The stokes axis is important, otherwise the cache is set up
752 : // (by the TiledStMan) such that it can hold only 1 stokes
753 : // with the result that iterating is tremendously slow.
754 : // We also need to cast the const away from _image.
755 30 : const IPosition mainShape = _image->shape();
756 15 : uInt nrtiles = (1 + (mainShape(fAxis)-1) / tileShape(fAxis)) *
757 15 : (1 + (mainShape(sAxis)-1) / tileShape(sAxis));
758 15 : auto* mainImagePtr = const_cast<ImageInterface<Float>*>(_image.get());
759 15 : mainImagePtr->setCacheSizeInTiles (nrtiles);
760 30 : String posString;
761 15 : auto ok = false;
762 30 : IPosition shp;
763 6015 : for (it.reset(); ! it.atEnd(); it++) {
764 : // Find rotation measure for this line
765 6000 : ok = _findRotationMeasure(
766 : rm, rmErr, pa0, pa0Err, rChiSq, nTurns, sortidx, wsqsort,
767 12000 : it.vectorCursor(), it.getMask(false),
768 12000 : paerr.getSlice(it.position(),it.cursorShape()),
769 : rmFg, rmMax, maxPaErr, posString
770 : );
771 : // Plonk values into output image. This is slow and clunky, but
772 : // should be relatively fast c.f. the fitting. Could be reimplemented
773 : // with LatticeApply if need be. Buffering is hard because the
774 : // navigator doesn't take a regular path. If I used a LatticeStepper
775 : // instead, the path would be regular and then I could buffer, but then
776 : // the iteration would be less efficient !!!
777 6000 : j = k = l = m = 0;
778 30000 : for (Int i=0; i<Int(it.position().size()); ++i) {
779 24000 : if (doRM && i != fAxis && i != sAxis) {
780 12000 : whereRM(j) = it.position()[i];
781 12000 : ++j;
782 : }
783 24000 : if (doPA && i != fAxis) {
784 1200 : wherePA(k) = it.position()[i];
785 1200 : ++k;
786 : }
787 24000 : if (doNTurns && i != fAxis && i != sAxis) {
788 0 : whereNTurns[l] = it.position()[i];
789 0 : ++l;
790 : }
791 24000 : if (doChiSq && i != fAxis && i != sAxis) {
792 0 : whereChiSq[m] = it.position()[i];
793 0 : ++m;
794 : }
795 : }
796 6000 : if (isMaskedRM) {
797 6000 : tmpMaskRM.set(ok);
798 6000 : outRMMaskPtr->putSlice(tmpMaskRM, whereRM);
799 : }
800 6000 : if (isMaskedRMErr) {
801 0 : tmpMaskRM.set(ok);
802 0 : outRMErrMaskPtr->putSlice(tmpMaskRM, whereRM);
803 : }
804 6000 : if (isMaskedPa0) {
805 400 : tmpMaskPA.set(ok);
806 400 : outPa0MaskPtr->putSlice(tmpMaskPA, wherePA);
807 : }
808 6000 : if (isMaskedPa0Err) {
809 0 : tmpMaskPA.set(ok);
810 0 : outPa0ErrMaskPtr->putSlice(tmpMaskPA, wherePA);
811 : }
812 6000 : if (isMaskedNTurns) {
813 0 : tmpMaskNTurns.set(ok);
814 0 : outNTurnsMaskPtr->putSlice(tmpMaskNTurns, whereNTurns);
815 : }
816 6000 : if (isMaskedChiSq) {
817 0 : tmpMaskChiSq.set(ok);
818 0 : outChiSqMaskPtr->putSlice(tmpMaskChiSq, whereChiSq);
819 : }
820 : // If the output value is masked, the value itself is 0
821 6000 : if (rmOutPtr) {
822 6000 : tmpValueRM.set(rm);
823 6000 : rmOutPtr->putSlice(tmpValueRM, whereRM);
824 : }
825 6000 : if (rmOutErrorPtr) {
826 0 : tmpValueRM.set(rmErr);
827 0 : rmOutErrorPtr->putSlice(tmpValueRM, whereRM);
828 : }
829 : // Position angles in degrees
830 6000 : if (pa0OutPtr) {
831 400 : tmpValuePA.set(pa0*degPerRad);
832 400 : pa0OutPtr->putSlice(tmpValuePA, wherePA);
833 : }
834 6000 : if (pa0OutErrorPtr) {
835 0 : tmpValuePA.set(pa0Err*degPerRad);
836 0 : pa0OutErrorPtr->putSlice(tmpValuePA, wherePA);
837 : }
838 : // Number of turns and chi sq
839 6000 : if (nTurnsOutPtr) {
840 0 : tmpValueNTurns.set(nTurns);
841 0 : nTurnsOutPtr->putSlice(tmpValueNTurns, whereNTurns);
842 : }
843 6000 : if (chiSqOutPtr) {
844 0 : tmpValueChiSq.set(rChiSq);
845 0 : chiSqOutPtr->putSlice(tmpValueChiSq, whereChiSq);
846 : }
847 6000 : if (showProgress) {
848 6000 : pProgressMeter->update(Double(it.nsteps()));
849 : }
850 : }
851 : // Clear the cache of the main image again.
852 15 : mainImagePtr->clearCache();
853 15 : }
854 :
855 32 : IPosition ImagePolarimetry::rotationMeasureShape(
856 : CoordinateSystem& cSys, Int& fAxis, Int& sAxis, LogIO&, Int spectralAxis
857 : ) const {
858 64 : const auto cSys0 = coordinates();
859 : Int spectralCoord;
860 32 : _findFrequencyAxis(spectralCoord, fAxis, cSys0, spectralAxis);
861 32 : Int afterCoord = -1;
862 32 : auto stokesCoord = cSys0.findCoordinate(Coordinate::STOKES, afterCoord);
863 64 : auto pixelAxes = cSys0.pixelAxes(stokesCoord);
864 32 : sAxis = pixelAxes(0);
865 : // What shape should the image be ? Frequency and stokes axes should be gone.
866 64 : const auto shape0 = shape();
867 32 : IPosition shape(shape0.size()-2);
868 32 : Int j = 0;
869 160 : for (Int i=0; i<Int(shape0.size()); ++i) {
870 128 : if (i != fAxis && i != sAxis) {
871 64 : shape[j] = shape0[i];
872 64 : ++j;
873 : }
874 : }
875 64 : CoordinateSystem tmp;
876 32 : cSys = tmp;
877 128 : for (Int i=0; i<Int(cSys0.nCoordinates()); ++i) {
878 96 : if (i != spectralCoord && i != stokesCoord) {
879 32 : cSys.addCoordinate(cSys0.coordinate(i));
880 : }
881 : }
882 64 : return shape;
883 : }
884 :
885 32 : IPosition ImagePolarimetry::positionAngleShape(
886 : CoordinateSystem& cSys, Int& fAxis, Int& sAxis, LogIO&, Int spectralAxis
887 : ) const {
888 64 : CoordinateSystem cSys0 = coordinates();
889 32 : Int spectralCoord = -1;
890 32 : _findFrequencyAxis (spectralCoord, fAxis, cSys0, spectralAxis);
891 32 : Int afterCoord = -1;
892 32 : Int stokesCoord = cSys0.findCoordinate(Coordinate::STOKES, afterCoord);
893 64 : Vector<Int> pixelAxes = cSys0.pixelAxes(stokesCoord);
894 32 : sAxis = pixelAxes(0);
895 32 : _fiddleStokesCoordinate(cSys0, Stokes::Pangle);
896 64 : CoordinateSystem tmp;
897 32 : cSys = tmp;
898 128 : for (Int i=0; i<Int(cSys0.nCoordinates()); ++i) {
899 96 : if (i != spectralCoord) {
900 64 : cSys.addCoordinate(cSys0.coordinate(i));
901 : }
902 : }
903 : // What shape should the image be ? Frequency axis should be gone.
904 : // and Stokes length 1
905 64 : const auto shape0 = ImagePolarimetry::shape();
906 32 : IPosition shape(shape0.size()-1);
907 32 : Int j = 0;
908 160 : for (Int i=0; i<Int(shape0.size()); ++i) {
909 128 : if (i == sAxis) {
910 32 : shape[j] = 1;
911 32 : ++j;
912 : }
913 : else {
914 96 : if (i != fAxis) {
915 64 : shape[j] = shape0[i];
916 64 : ++j;
917 : }
918 : }
919 : }
920 64 : return shape;
921 : }
922 :
923 0 : ImageExpr<Float> ImagePolarimetry::stokesI() const {
924 : return _makeStokesExpr(
925 0 : _stokes[ImagePolarimetry::I], Stokes::I, "StokesI"
926 0 : );
927 : }
928 :
929 0 : Float ImagePolarimetry::sigmaStokesI(Float clip) {
930 0 : return _sigma(ImagePolarimetry::I, clip);
931 : }
932 :
933 0 : ImageExpr<Float> ImagePolarimetry::stokesQ() const {
934 0 : return _makeStokesExpr(_stokes[ImagePolarimetry::Q], Stokes::Q, "StokesQ");
935 : }
936 :
937 0 : Float ImagePolarimetry::sigmaStokesQ(Float clip) {
938 0 : return _sigma(ImagePolarimetry::Q, clip);
939 : }
940 :
941 0 : ImageExpr<Float> ImagePolarimetry::stokesU() const {
942 0 : return _makeStokesExpr(_stokes[ImagePolarimetry::U], Stokes::U, "StokesU");
943 : }
944 :
945 0 : Float ImagePolarimetry::sigmaStokesU(Float clip) {
946 0 : return _sigma(ImagePolarimetry::U, clip);
947 : }
948 :
949 0 : ImageExpr<Float> ImagePolarimetry::stokesV() const {
950 0 : return _makeStokesExpr(_stokes[ImagePolarimetry::V], Stokes::V, "StokesV");
951 : }
952 :
953 0 : Float ImagePolarimetry::sigmaStokesV(Float clip) {
954 0 : return _sigma(ImagePolarimetry::V, clip);
955 : }
956 :
957 0 : ImageExpr<Float> ImagePolarimetry::stokes(
958 : ImagePolarimetry::StokesTypes stokes
959 : ) const {
960 0 : const auto type = _stokesType(stokes);
961 0 : return _makeStokesExpr(_stokes[stokes], type, _stokesName(stokes));
962 : }
963 :
964 0 : Float ImagePolarimetry::sigmaStokes(
965 : ImagePolarimetry::StokesTypes stokes, Float clip
966 : ) {
967 0 : return _sigma(stokes, clip);
968 : }
969 :
970 0 : void ImagePolarimetry::summary(LogIO& os) const {
971 0 : ImageSummary<Float> s(*_image);
972 0 : s.list(os);
973 0 : }
974 :
975 0 : ImageExpr<Float> ImagePolarimetry::totPolInt(
976 : Bool debias, Float clip, Float sigma
977 : ) {
978 0 : LogIO os(LogOrigin("ImagePolarimetry", __func__, WHERE));
979 : Bool doLin, doCirc;
980 0 : _setDoLinDoCirc(doLin, doCirc);
981 0 : auto node = _makePolIntNode(os, debias, clip, sigma, doLin, doCirc);
982 0 : LatticeExpr<Float> le(node);
983 0 : ImageExpr<Float> ie(le, String("totalPolarizedIntensity"));
984 : // Dodgy. The beam is now rectified
985 0 : ie.setUnits(_image->units());
986 0 : StokesTypes stokes = _stokes[Q] ? Q : _stokes[U] ? U : V;
987 0 : _setInfo(ie, stokes);
988 0 : _fiddleStokesCoordinate(ie, Stokes::Ptotal);
989 0 : return ie;
990 : }
991 :
992 0 : Float ImagePolarimetry::sigmaTotPolInt(Float clip, Float sigma) {
993 : // sigma_P = sigma_QUV
994 0 : LogIO os(LogOrigin("ImagePolarimetry", __func__, WHERE));
995 : Bool doLin, doCirc;
996 0 : _setDoLinDoCirc(doLin, doCirc);
997 0 : Float sigma2 = sigma > 0 ? sigma : this->sigma(clip);
998 0 : return sigma2;
999 : }
1000 :
1001 :
1002 0 : IPosition ImagePolarimetry::singleStokesShape(
1003 : CoordinateSystem& cSys, Stokes::StokesTypes type
1004 : ) const {
1005 : // We know the image has a Stokes coordinate or it
1006 : // would have failed at construction
1007 0 : auto cSys0 = _image->coordinates();
1008 0 : _fiddleStokesCoordinate(cSys0, type);
1009 0 : cSys = cSys0;
1010 0 : Int afterCoord = -1;
1011 0 : const auto iStokes = cSys0.findCoordinate(Coordinate::STOKES, afterCoord);
1012 0 : const auto pixelAxes = cSys0.pixelAxes(iStokes);
1013 0 : auto shape = _image->shape();
1014 0 : shape[pixelAxes[0]] = 1;
1015 0 : return shape;
1016 : }
1017 :
1018 0 : ImageExpr<Float> ImagePolarimetry::depolarizationRatio(
1019 : const ImageInterface<Float>& im1, const ImageInterface<Float>& im2,
1020 : Bool debias, Float clip, Float sigma
1021 : ) {
1022 0 : ImagePolarimetry p1(im1);
1023 0 : ImagePolarimetry p2(im2);
1024 0 : ImageExpr<Float> m1(p1.fracLinPol(debias, clip, sigma));
1025 0 : ImageExpr<Float> m2(p2.fracLinPol(debias, clip, sigma));
1026 0 : LatticeExprNode n1(m1/m2);
1027 0 : LatticeExpr<Float> le(n1);
1028 0 : ImageExpr<Float> depol(le, "DepolarizationRatio");
1029 0 : return depol;
1030 : }
1031 :
1032 0 : ImageExpr<Float> ImagePolarimetry::sigmaDepolarizationRatio(
1033 : const ImageInterface<Float>& im1, const ImageInterface<Float>& im2,
1034 : Bool debias, Float clip, Float sigma
1035 : ) {
1036 0 : Vector<StokesTypes> types(3);
1037 0 : types[0] = I;
1038 0 : types[1] = Q;
1039 0 : types[2] = U;
1040 0 : _checkBeams(im1, im2, types);
1041 0 : ImagePolarimetry p1(im1);
1042 0 : ImagePolarimetry p2(im2);
1043 0 : ImageExpr<Float> m1 = p1.fracLinPol(debias, clip, sigma);
1044 0 : ImageExpr<Float> sm1 = p1.sigmaFracLinPol(clip, sigma);
1045 0 : ImageExpr<Float> m2 = p2.fracLinPol(debias, clip, sigma);
1046 0 : ImageExpr<Float> sm2 = p2.sigmaFracLinPol(clip, sigma);
1047 0 : LatticeExprNode n0(m1/m2);
1048 0 : LatticeExprNode n1(sm1*sm1/m1/m1);
1049 0 : LatticeExprNode n2(sm2*sm2/m2/m2);
1050 0 : LatticeExprNode n3(n0 * sqrt(n1+n2));
1051 0 : LatticeExpr<Float> le(n3);
1052 0 : ImageExpr<Float> sigmaDepol(le, "DepolarizationRatioError");
1053 0 : return sigmaDepol;
1054 : }
1055 :
1056 44 : void ImagePolarimetry::_cleanup() {
1057 44 : _image.reset();
1058 220 : for (uInt i=0; i<4; ++i) {
1059 176 : delete _stokes[i];
1060 176 : _stokes[i] = nullptr;
1061 176 : delete _stokesStats[i];
1062 176 : _stokesStats[i] = nullptr;
1063 : }
1064 44 : delete _fitter;
1065 44 : _fitter = nullptr;
1066 44 : }
1067 :
1068 64 : void ImagePolarimetry::_findFrequencyAxis(
1069 : Int& spectralCoord, Int& fAxis, const CoordinateSystem& cSys,
1070 : Int spectralAxis
1071 : ) const {
1072 192 : LogIO os(LogOrigin("ImagePolarimetry", __func__, WHERE));
1073 64 : spectralCoord = -1;
1074 64 : fAxis = -1;
1075 64 : if (spectralAxis >= 0) {
1076 0 : ThrowIf(
1077 : spectralAxis >= Int(cSys.nPixelAxes()),
1078 : "Illegal spectral axis " + String::toString(spectralAxis) +" given"
1079 : );
1080 0 : fAxis = spectralAxis;
1081 : Int axisInCoordinate;
1082 0 : cSys.findPixelAxis(spectralCoord, axisInCoordinate, fAxis);
1083 : // Check coordinate type is one of expected types
1084 0 : ThrowIf(
1085 : ! (
1086 : cSys.type(spectralCoord)==Coordinate::TABULAR
1087 : || cSys.type(spectralCoord)==Coordinate::LINEAR
1088 : || cSys.type(spectralCoord)==Coordinate::SPECTRAL
1089 : ),
1090 : "The specified axis of type " + cSys.showType(spectralCoord)
1091 : + " cannot be a frequency axis"
1092 : );
1093 : }
1094 : else {
1095 64 : spectralCoord = _findSpectralCoordinate(cSys, os, false);
1096 64 : if (spectralCoord < 0) {
1097 0 : for (uInt i=0; i<cSys.nCoordinates(); ++i) {
1098 0 : if (
1099 0 : cSys.type(i)==Coordinate::TABULAR
1100 0 : || cSys.type(i)==Coordinate::LINEAR
1101 : ) {
1102 0 : const auto axisNames = cSys.coordinate(i).worldAxisNames();
1103 0 : String tmp = axisNames[0];
1104 0 : tmp.upcase();
1105 0 : if (tmp.contains(String("FREQ"))) {
1106 0 : spectralCoord = i;
1107 0 : break;
1108 : }
1109 : }
1110 : }
1111 : }
1112 64 : ThrowIf(
1113 : spectralCoord < 0,
1114 : "Cannot find SpectralCoordinate in this image"
1115 : );
1116 64 : fAxis = cSys.pixelAxes(spectralCoord)[0];
1117 : }
1118 64 : }
1119 :
1120 44 : void ImagePolarimetry::_findStokes() {
1121 132 : LogIO os(LogOrigin("ImagePolarimetry", __func__, WHERE));
1122 44 : const CoordinateSystem& cSys = _image->coordinates();
1123 44 : Int polAxisNum = cSys.polarizationAxisNumber();
1124 44 : ThrowIf(
1125 : polAxisNum < 0, "There is no Stokes Coordinate in this image"
1126 : );
1127 44 : const uInt ndim = _image->ndim();
1128 88 : auto shape = _image->shape();
1129 88 : IPosition blc(ndim,0);
1130 88 : auto trc = shape - 1;
1131 220 : for (const auto& kv : polMap) {
1132 176 : const auto pix = cSys.stokesPixelNumber(kv.second);
1133 176 : if (pix >= 0) {
1134 136 : _stokes[kv.first] = _makeSubImage(blc, trc, polAxisNum, pix);
1135 : }
1136 : }
1137 44 : if((_stokes[Q] && ! _stokes[U]) || (! _stokes[Q] && _stokes[U])) {
1138 5 : _cleanup();
1139 5 : ThrowCc(
1140 : "This Stokes coordinate has only one of Q and U. This is not useful"
1141 : );
1142 : }
1143 39 : if (! (_stokes[Q] || _stokes[U] || _stokes[V])) {
1144 1 : _cleanup();
1145 1 : ThrowCc("This image has no Stokes Q, U, nor V. This is not useful");
1146 : }
1147 38 : }
1148 :
1149 32 : void ImagePolarimetry::_fiddleStokesCoordinate(
1150 : ImageInterface<Float>& im, Stokes::StokesTypes type
1151 : ) const {
1152 64 : CoordinateSystem cSys = im.coordinates();
1153 32 : _fiddleStokesCoordinate(cSys, type);
1154 32 : im.setCoordinateInfo(cSys);
1155 32 : }
1156 :
1157 64 : void ImagePolarimetry::_fiddleStokesCoordinate(
1158 : CoordinateSystem& cSys, Stokes::StokesTypes type
1159 : ) const {
1160 64 : Int afterCoord = -1;
1161 64 : Int iStokes = cSys.findCoordinate(Coordinate::STOKES, afterCoord);
1162 128 : const Vector<Int> which(1, Int(type));
1163 128 : StokesCoordinate stokes(which);
1164 64 : cSys.replaceCoordinate(stokes, iStokes);
1165 64 : }
1166 :
1167 0 : void ImagePolarimetry::_fiddleStokesCoordinate(
1168 : ImageInterface<Complex>& ie, Stokes::StokesTypes type
1169 : ) const {
1170 0 : CoordinateSystem cSys = ie.coordinates();
1171 0 : _fiddleStokesCoordinate(cSys, type);
1172 0 : ie.setCoordinateInfo(cSys);
1173 0 : }
1174 :
1175 0 : void ImagePolarimetry::_fiddleTimeCoordinate(
1176 : ImageInterface<Complex>& ie, const Quantum<Double>& f, Int coord
1177 : ) const {
1178 0 : LogIO os(LogOrigin("ImagePolarimetry", __func__, WHERE));
1179 0 : CoordinateSystem cSys = ie.coordinates();
1180 0 : unique_ptr<Coordinate> pC(cSys.coordinate(coord).clone());
1181 0 : AlwaysAssert(pC->nPixelAxes()==1,AipsError);
1182 0 : AlwaysAssert(pC->type()==Coordinate::LINEAR,AipsError);
1183 0 : auto axisUnits = pC->worldAxisUnits();
1184 0 : axisUnits = String("s");
1185 0 : ThrowIf(
1186 : ! pC->setWorldAxisUnits(axisUnits),
1187 : "Failed to set TimeCoordinate units to seconds because "
1188 : + pC->errorMessage()
1189 : );
1190 : // Find factor to convert from time (s) to rad/m/m
1191 0 : auto inc = pC->increment();
1192 0 : const auto ff = f.getValue(Unit("Hz"));
1193 0 : const auto lambda = QC::c( ).getValue(Unit("m/s")) / ff;
1194 0 : const auto fac = -C::pi * ff / 2.0 / lambda / lambda;
1195 0 : inc *= fac;
1196 0 : Vector<String> axisNames(1);
1197 0 : axisNames = String("RotationMeasure");
1198 0 : axisUnits = String("rad/m/m");
1199 0 : Vector<Double> refVal(1,0.0);
1200 : LinearCoordinate lC(
1201 : axisNames, axisUnits, refVal, inc,
1202 0 : pC->linearTransform().copy(), pC->referencePixel().copy()
1203 0 : );
1204 0 : cSys.replaceCoordinate(lC, coord);
1205 0 : ie.setCoordinateInfo(cSys);
1206 0 : }
1207 :
1208 0 : Quantum<Double> ImagePolarimetry::_findCentralFrequency(
1209 : const Coordinate& coord, Int shape
1210 : ) const {
1211 0 : AlwaysAssert(coord.nPixelAxes()==1,AipsError);
1212 0 : Vector<Double> pixel(1);
1213 0 : Vector<Double> world;
1214 0 : pixel(0) = Double(shape - 1) / 2.0;
1215 0 : ThrowIf(
1216 : ! coord.toWorld(world, pixel),
1217 : "Failed to convert pixel to world for SpectralCoordinate because "
1218 : + coord.errorMessage()
1219 : );
1220 0 : const auto units = coord.worldAxisUnits();
1221 0 : return Quantum<Double>(world(0), units(0));
1222 : }
1223 :
1224 64 : Int ImagePolarimetry::_findSpectralCoordinate(
1225 : const CoordinateSystem& cSys, LogIO& os, Bool fail
1226 : ) const {
1227 64 : Int afterCoord = -1;
1228 64 : Int coord = cSys.findCoordinate(Coordinate::SPECTRAL, afterCoord);
1229 64 : ThrowIf(coord < 0 && fail, "No spectral coordinate in this image");
1230 64 : if (afterCoord>0) {
1231 : os << LogIO::WARN << "This image has more than one spectral "
1232 0 : << "coordinate; only first used" << LogIO::POST;
1233 : }
1234 64 : return coord;
1235 : }
1236 :
1237 6000 : Bool ImagePolarimetry::_findRotationMeasure(
1238 : Float& rmFitted, Float& rmErrFitted, Float& pa0Fitted, Float& pa0ErrFitted,
1239 : Float& rChiSqFitted, Float& nTurns, const Vector<uInt>& sortidx,
1240 : const Vector<Float>& wsq2, const Vector<Float>& pa2,
1241 : const Array<Bool>& paMask2, const Array<Float>& paerr2, Float rmFg,
1242 : Float rmMax, Float maxPaErr, const String& posString
1243 : ) {
1244 : // wsq is lambda squared in m**2 in increasing wavelength order
1245 : // pa is position angle in radians
1246 : // paerr is pa error in radians
1247 : // maxPaErr is maximum tolerated error in position angle
1248 : // rmfg is a user specified foreground RM rad/m/m
1249 : // rmmax is a user specified maximum RM
1250 6000 : static Vector<Float> paerr;
1251 6000 : static Vector<Float> pa;
1252 6000 : static Vector<Float> wsq;
1253 : // Abandon if less than 2 points
1254 6000 : uInt n = sortidx.size();
1255 6000 : if (n<2) {
1256 0 : return false;
1257 : }
1258 6000 : rmFitted = rmErrFitted = pa0Fitted = pa0ErrFitted = rChiSqFitted = 0.0;
1259 : // Sort into decreasing frequency order and correct for foreground rotation
1260 : // Remember wsq already sorted. Discard points that are too noisy or masked
1261 12000 : const Vector<Float>& paerr1(paerr2.nonDegenerate(0));
1262 12000 : const Vector<Bool>& paMask1(paMask2.nonDegenerate(0));
1263 6000 : paerr.resize(n);
1264 6000 : pa.resize(n);
1265 6000 : wsq.resize(n);
1266 6000 : uInt j = 0;
1267 144000 : for (uInt i=0; i<n; ++i) {
1268 138000 : if (abs(paerr1(sortidx(i))) < maxPaErr && paMask1(sortidx(i))) {
1269 130000 : pa(j) = pa2(sortidx(i)) - rmFg*wsq2(i);
1270 130000 : paerr(j) = paerr1(sortidx(i));
1271 130000 : wsq(j) = wsq2(i);
1272 130000 : ++j;
1273 : }
1274 : }
1275 6000 : n = j;
1276 6000 : if (n<=1) {
1277 400 : return false;
1278 : }
1279 5600 : pa.resize(n, true);
1280 5600 : paerr.resize(n, true);
1281 5600 : wsq.resize(n, true);
1282 : // Treat supplementary and primary points separately
1283 : Bool ok = n == 2
1284 5600 : ? _rmSupplementaryFit(
1285 : nTurns, rmFitted, rmErrFitted, pa0Fitted, pa0ErrFitted,
1286 : rChiSqFitted, wsq, pa, paerr
1287 : )
1288 5600 : : _rmPrimaryFit(
1289 : nTurns, rmFitted, rmErrFitted, pa0Fitted, pa0ErrFitted,
1290 : rChiSqFitted, wsq, pa, paerr, rmMax, /*plotter,*/ posString
1291 5600 : );
1292 : // Put position angle into the range 0->pi
1293 5600 : static MVAngle tmpMVA1;
1294 5600 : if (ok) {
1295 5600 : MVAngle tmpMVA0(pa0Fitted);
1296 5600 : tmpMVA1 = tmpMVA0.binorm(0.0);
1297 5600 : pa0Fitted = tmpMVA1.radian();
1298 : // Add foreground back on
1299 5600 : rmFitted += rmFg;
1300 : }
1301 5600 : return ok;
1302 : }
1303 :
1304 16 : void ImagePolarimetry::_hasQU () const {
1305 16 : ThrowIf(
1306 : ! (_stokes[Q] && _stokes[U]),
1307 : "This image does not have Stokes Q and U which are "
1308 : "required for this function"
1309 : );
1310 16 : }
1311 :
1312 0 : ImageExpr<Float> ImagePolarimetry::_makeStokesExpr(
1313 : ImageInterface<Float>* imPtr, Stokes::StokesTypes type, const String& name
1314 : ) const {
1315 0 : ThrowIf(! imPtr, "This image does not have Stokes " + Stokes::name(type));
1316 0 : LatticeExprNode node(*imPtr);
1317 0 : LatticeExpr<Float> le(node);
1318 0 : ImageExpr<Float> ie(le, name);
1319 0 : ie.setUnits(_image->units());
1320 0 : _fiddleStokesCoordinate(ie, type);
1321 0 : return ie;
1322 : }
1323 :
1324 136 : ImageInterface<Float>* ImagePolarimetry::_makeSubImage(
1325 : IPosition& blc, IPosition& trc, Int axis, Int pix
1326 : ) const {
1327 136 : blc[axis] = pix;
1328 136 : trc[axis] = pix;
1329 272 : LCSlicer slicer(blc, trc, RegionType::Abs);
1330 136 : ImageRegion region(slicer);
1331 272 : return new SubImage<Float>(*_image, region);
1332 : }
1333 :
1334 0 : LatticeExprNode ImagePolarimetry::_makePolIntNode(
1335 : LogIO& os, Bool debias, Float clip, Float sigma, Bool doLin, Bool doCirc
1336 : ) {
1337 0 : LatticeExprNode linNode, circNode, node;
1338 0 : Float sigma2 = debias ? (sigma > 0.0 ? sigma : this->sigma(clip)) : 0.0;
1339 0 : if (doLin) {
1340 0 : linNode = LatticeExprNode(pow(*_stokes[U], 2) + pow(*_stokes[Q], 2));
1341 : }
1342 0 : if (doCirc) {
1343 0 : circNode = LatticeExprNode(pow(*_stokes[V], 2));
1344 : }
1345 0 : Float sigmasq = sigma2 * sigma2;
1346 0 : if (doLin && doCirc) {
1347 0 : node = linNode + circNode;
1348 0 : if (debias) {
1349 0 : node = node - LatticeExprNode(sigmasq);
1350 0 : os << LogIO::NORMAL << "Debiasing with sigma = " << sqrt(sigmasq)
1351 0 : << LogIO::POST;
1352 : }
1353 : }
1354 0 : else if (doLin) {
1355 0 : node = linNode;
1356 0 : if (debias) {
1357 0 : node = node - LatticeExprNode(sigmasq);
1358 0 : os << LogIO::NORMAL << "Debiasing with sigma = " << sqrt(sigmasq)
1359 0 : << LogIO::POST;
1360 : }
1361 : }
1362 0 : else if (doCirc) {
1363 0 : node = circNode;
1364 0 : if (debias) {
1365 0 : node = node - LatticeExprNode(sigmasq);
1366 0 : os << LogIO::NORMAL << "Debiasing with sigma = " << sqrt(sigmasq)
1367 0 : << LogIO::POST;
1368 : }
1369 : }
1370 0 : return LatticeExprNode(sqrt(node));
1371 : }
1372 :
1373 5600 : Bool ImagePolarimetry::_rmPrimaryFit(
1374 : Float& nTurns, Float& rmFitted, Float& rmErrFitted, Float& pa0Fitted,
1375 : Float& pa0ErrFitted, Float& rChiSqFitted, const Vector<Float>& wsq,
1376 : const Vector<Float>& pa, const Vector<Float>& paerr, Float rmMax,
1377 : const String&
1378 : ) {
1379 5600 : static Vector<Float> plotPA;
1380 5600 : static Vector<Float> plotPAErr;
1381 5600 : static Vector<Float> plotPAErrY1;
1382 5600 : static Vector<Float> plotPAErrY2;
1383 5600 : static Vector<Float> plotPAFit;
1384 : // Assign position angle to longest wavelength consistent with RM < RMMax
1385 5600 : const uInt n = wsq.size();
1386 5600 : Double dwsq = wsq(n-1) - wsq(0);
1387 5600 : Float ppa = abs(rmMax)*dwsq + pa(0);
1388 5600 : Float diff = ppa - pa(n-1);
1389 5600 : Float t = diff >= 0 ? 0.5 : -0.5;
1390 5600 : Int maxnpi = Int(diff/C::pi + t);
1391 5600 : ppa = -abs(rmMax)*dwsq + pa(0);
1392 5600 : diff = ppa - pa(n-1);
1393 5600 : t = diff >= 0 ? 0.5 : -0.5;
1394 5600 : Int minnpi = Int(diff/C::pi + t);
1395 : // Loop over range of n*pi ambiguity
1396 11200 : Vector<Float> fitpa(n);
1397 11200 : Vector<Float> pars;
1398 5600 : Float chiSq = 1e30;
1399 41936 : for (Int h=minnpi; h<=maxnpi; ++h) {
1400 36336 : fitpa[n-1] = pa[n-1] + C::pi*h;
1401 36336 : Float rm0 = (fitpa(n-1) - pa(0))/dwsq;
1402 : // Assign position angles to remaining wavelengths
1403 708384 : for (uInt k=1; k<n-1; ++k) {
1404 672048 : ppa = pa[0] + rm0*(wsq[k]-wsq[0]);
1405 672048 : diff = ppa - pa[k];
1406 672048 : t = diff >= 0 ? 0.5 : -0.5;
1407 672048 : Int npi = Int(diff/C::pi + t);
1408 672048 : fitpa[k] = pa[k] + npi*C::pi;
1409 : }
1410 36336 : fitpa[0] = pa[0];
1411 : // Do least squares fit
1412 36336 : if (!_rmLsqFit (pars, wsq, fitpa, paerr)) {
1413 0 : return false;
1414 : }
1415 36336 : if (pars[4] < chiSq) {
1416 7199 : chiSq = pars[4];
1417 7199 : nTurns = h; // Number of turns
1418 7199 : rmFitted = pars[0]; // Fitted RM
1419 7199 : rmErrFitted = pars[1]; // Error in RM
1420 7199 : pa0Fitted = pars[2]; // Fitted intrinsic angle
1421 7199 : pa0ErrFitted = pars[3]; // Error in angle
1422 7199 : rChiSqFitted = pars[4]; // Recued chi squared
1423 7199 : if (n > 2) {
1424 7199 : rChiSqFitted /= Float(n - 2);
1425 : }
1426 : }
1427 : }
1428 5600 : return true;
1429 : }
1430 :
1431 0 : Bool ImagePolarimetry::_rmSupplementaryFit(
1432 : Float& nTurns, Float& rmFitted, Float& rmErrFitted, Float& pa0Fitted,
1433 : Float& pa0ErrFitted, Float& rChiSqFitted, const Vector<Float>& wsq,
1434 : const Vector<Float>& pa, const Vector<Float>& paerr
1435 : ) {
1436 : // For supplementary points find lowest residual RM
1437 0 : const uInt n = wsq.size();
1438 0 : Float absRM = 1e30;
1439 0 : Vector<Float> fitpa(pa.copy());
1440 0 : Vector<Float> pars;
1441 0 : for (Int i=-2; i<3; ++i) {
1442 0 : fitpa[n-1] = pa[n-1] + C::pi*i;
1443 : // Do least squares fit
1444 0 : if (! _rmLsqFit (pars, wsq, fitpa, paerr)) {
1445 0 : return false;
1446 : }
1447 : // Save solution with lowest absolute RM
1448 0 : if (abs(pars[0]) < absRM) {
1449 0 : absRM = abs(pars[0]);
1450 0 : nTurns = i; // nTurns
1451 0 : rmFitted = pars(0); // Fitted RM
1452 0 : rmErrFitted = pars[1]; // Error in RM
1453 0 : pa0Fitted = pars[2]; // Fitted intrinsic angle
1454 0 : pa0ErrFitted = pars[3]; // Error in angle
1455 0 : rChiSqFitted = pars[4]; // Reduced chi squared
1456 0 : if (n > 2) {
1457 0 : rChiSqFitted /= Float(n - 2);
1458 : }
1459 : }
1460 : }
1461 0 : return true;
1462 : }
1463 :
1464 36336 : Bool ImagePolarimetry::_rmLsqFit(
1465 : Vector<Float>& pars, const Vector<Float>& wsq, const Vector<Float> pa,
1466 : const Vector<Float>& paerr
1467 : ) const {
1468 : // Perform fit on unmasked data
1469 36336 : static Vector<Float> solution;
1470 : try {
1471 36336 : solution = _fitter->fit(wsq, pa, paerr);
1472 0 : } catch (AipsError x) {
1473 0 : return false;
1474 : }
1475 36336 : const Vector<Double>& cv = _fitter->compuCovariance().diagonal();
1476 36336 : pars.resize(5);
1477 36336 : pars[0] = solution[1];
1478 36336 : pars[1] = sqrt(cv[1]);
1479 36336 : pars[2] = solution[0];
1480 36336 : pars[3] = sqrt(cv[0]);
1481 36336 : pars[4] = _fitter->chiSquare();
1482 36336 : return true;
1483 : }
1484 :
1485 0 : String ImagePolarimetry::_stokesName(
1486 : ImagePolarimetry::StokesTypes index
1487 : ) const {
1488 0 : const auto iter = polMap.find(index);
1489 0 : return (iter == polMap.end()) ? "??" : iter->second;
1490 : }
1491 :
1492 0 : Stokes::StokesTypes ImagePolarimetry::_stokesType (ImagePolarimetry::StokesTypes index) const
1493 : {
1494 0 : if (index==ImagePolarimetry::I) {
1495 0 : return Stokes::I;
1496 0 : } else if (index==ImagePolarimetry::Q) {
1497 0 : return Stokes::Q;
1498 0 : } else if (index==ImagePolarimetry::U) {
1499 0 : return Stokes::U;
1500 0 : } else if (index==ImagePolarimetry::V) {
1501 0 : return Stokes::V;
1502 : } else {
1503 0 : return Stokes::Undefined;
1504 : }
1505 : }
1506 :
1507 :
1508 12 : Float ImagePolarimetry::_sigma(
1509 : ImagePolarimetry::StokesTypes index, Float clip
1510 : ) {
1511 12 : Float clip2 = clip == 0 ? 10.0 : abs(clip);
1512 12 : if (clip2 != _oldClip && _stokesStats[index]) {
1513 0 : delete _stokesStats[index];
1514 0 : _stokesStats[index] = 0;
1515 : }
1516 12 : if (! _stokesStats[index]) {
1517 : // Find sigma for all points inside +/- clip-sigma of the mean
1518 : // More joys of LEL
1519 12 : const ImageInterface<Float>* p = _stokes[index];
1520 24 : LatticeExprNode n1 (*p);
1521 36 : LatticeExprNode n2 (n1[abs(n1-mean(n1)) < clip2*stddev(n1)]);
1522 12 : LatticeExpr<Float> le(n2);
1523 12 : _stokesStats[index] = new LatticeStatistics<Float>(le, false, false);
1524 : }
1525 12 : Array<Float> sigmaA;
1526 12 : _stokesStats[index]->getConvertedStatistic(sigmaA, LatticeStatsBase::SIGMA);
1527 12 : ThrowIf(
1528 : sigmaA.empty(),
1529 : "No good points in clipped determination of the noise for the Stokes "
1530 : + _stokesName(index) + " image"
1531 : );
1532 12 : _oldClip = clip2;
1533 24 : return sigmaA(IPosition(1,0));
1534 : }
1535 :
1536 0 : void ImagePolarimetry::_subtractProfileMean(
1537 : ImageInterface<Float>& im, uInt axis
1538 : ) const {
1539 0 : const IPosition tileShape = im.niceCursorShape();
1540 0 : TiledLineStepper ts(im.shape(), tileShape, axis);
1541 0 : LatticeIterator<Float> it(im, ts);
1542 : Float dMean;
1543 0 : if (im.isMasked()) {
1544 0 : const Lattice<Bool>& mask = im.pixelMask();
1545 0 : for (it.reset(); !it.atEnd(); it++) {
1546 0 : const auto& p = it.cursor();
1547 0 : const auto& m = mask.getSlice(it.position(), it.cursorShape());
1548 0 : const MaskedArray<Float> ma(p, m, true);
1549 0 : dMean = mean(ma);
1550 0 : it.rwVectorCursor() -= dMean;
1551 : }
1552 : }
1553 : else {
1554 0 : for (it.reset(); !it.atEnd(); it++) {
1555 0 : dMean = mean(it.vectorCursor());
1556 0 : it.rwVectorCursor() -= dMean;
1557 : }
1558 : }
1559 0 : }
1560 :
1561 16 : Bool ImagePolarimetry::_dealWithMask(
1562 : Lattice<Bool>*& pMask, ImageInterface<Float>*& pIm, LogIO& os,
1563 : const String& type
1564 : ) const {
1565 16 : auto isMasked = pIm->isMasked();
1566 16 : if (! isMasked) {
1567 0 : if (pIm->canDefineRegion()) {
1568 0 : pIm->makeMask("mask0", true, true, true, true);
1569 0 : isMasked = true;
1570 : }
1571 : else {
1572 : os << LogIO::WARN << "Could not create a mask for the output "
1573 0 : << type << " image" << LogIO::POST;
1574 : }
1575 : }
1576 16 : if (isMasked) {
1577 16 : pMask = &(pIm->pixelMask());
1578 16 : if (! pMask->isWritable()) {
1579 0 : isMasked = false;
1580 : os << LogIO::WARN << "The output " << type
1581 0 : << " image has a mask but it's not writable" << LogIO::POST;
1582 : }
1583 : }
1584 16 : return isMasked;
1585 : }
1586 :
1587 38 : void ImagePolarimetry::_createBeamsEqMat() {
1588 38 : _beamsEqMat.assign(Matrix<Bool>(4, 4, false));
1589 38 : Bool hasMultiBeams = _image->imageInfo().hasMultipleBeams();
1590 190 : for (uInt i=0; i<4; ++i) {
1591 532 : for (uInt j=i; j<4; ++j) {
1592 380 : if (! (_stokes[i] && _stokes[j])) {
1593 94 : _beamsEqMat(i, j) = false;
1594 : }
1595 286 : else if (i == j) {
1596 126 : _beamsEqMat(i, j) = true;
1597 : }
1598 160 : else if (hasMultiBeams) {
1599 0 : _beamsEqMat(i, j) = (
1600 0 : _stokes[i]->imageInfo().getBeamSet()
1601 0 : == _stokes[j]->imageInfo().getBeamSet()
1602 : );
1603 : }
1604 : else {
1605 160 : _beamsEqMat(i, j) = true;
1606 : }
1607 380 : _beamsEqMat(j, i) = _beamsEqMat(i, j);
1608 :
1609 : }
1610 : }
1611 38 : }
1612 :
1613 48 : Bool ImagePolarimetry::_checkBeams(
1614 : const vector<StokesTypes>& stokes, Bool requireChannelEquality, Bool throws
1615 : ) const {
1616 96 : for (
1617 240 : auto iter = stokes.cbegin(); iter != stokes.cend(); ++iter
1618 : ) {
1619 144 : for (
1620 384 : auto jiter=iter; jiter!=stokes.cend(); ++jiter
1621 : ) {
1622 144 : if (iter == jiter) {
1623 96 : continue;
1624 : }
1625 48 : if (! _beamsEqMat(*iter, *jiter)) {
1626 0 : ThrowIf(
1627 : throws,
1628 : "Input image has multiple beams and the corresponding "
1629 : "beams for the stokes planes necessary for this "
1630 : "computation are not equal."
1631 : );
1632 0 : return False;
1633 : }
1634 : }
1635 : }
1636 48 : if (
1637 : requireChannelEquality
1638 0 : && _stokes[stokes[0]]->coordinates().hasSpectralAxis()
1639 48 : && _stokes[stokes[0]]->imageInfo().hasMultipleBeams()
1640 : ) {
1641 : const auto& beamSet
1642 0 : = _stokes[stokes[0]]->imageInfo().getBeamSet().getBeams();
1643 0 : auto start = beamSet.cbegin();
1644 0 : ++start;
1645 0 : for (auto iter=start; iter!=beamSet.cend(); ++iter) {
1646 0 : if (*iter != *(beamSet.begin())) {
1647 0 : ThrowIf(
1648 : throws,
1649 : "At least one beam in this image is not equal to all the "
1650 : "others along its spectral axis so this computation cannot "
1651 : "be performed"
1652 : );
1653 0 : return False;
1654 : }
1655 : }
1656 : }
1657 48 : return true;
1658 : }
1659 :
1660 0 : Bool ImagePolarimetry::_checkIQUBeams(
1661 : Bool requireChannelEquality, Bool throws
1662 : ) const {
1663 0 : static const vector<StokesTypes> types = {I, Q, U};
1664 0 : return _checkBeams(types, requireChannelEquality, throws);
1665 : }
1666 :
1667 0 : Bool ImagePolarimetry::_checkIVBeams(
1668 : Bool requireChannelEquality, Bool throws
1669 : ) const {
1670 0 : static const vector<StokesTypes> types = {I, V};
1671 0 : return _checkBeams(types, requireChannelEquality, throws);
1672 : }
1673 :
1674 48 : Bool ImagePolarimetry::_checkQUBeams(
1675 : Bool requireChannelEquality, Bool throws
1676 : ) const {
1677 48 : static const vector<StokesTypes> types = {Q, U};
1678 48 : return _checkBeams(types, requireChannelEquality, throws);
1679 : }
1680 :
1681 0 : void ImagePolarimetry::_checkBeams(
1682 : const ImagePolarimetry& im1, const ImagePolarimetry& im2,
1683 : const Vector<ImagePolarimetry::StokesTypes>& stokes
1684 : ) {
1685 0 : for (auto iter=stokes.cbegin(); iter!=stokes.cend(); iter++) {
1686 0 : ThrowIf(
1687 : im1.stokes(*iter).imageInfo().getBeamSet()
1688 : != im2.stokes(*iter).imageInfo().getBeamSet(),
1689 : "Beams or beamsets are not equal between the two images, caonnot "
1690 : "perform calculation"
1691 : );
1692 : }
1693 0 : }
1694 :
1695 : }
|
