Line data Source code
1 : //# ImagePolProxy.cc: C++ casa namespace proxy for ImagePol tool
2 : //# Copyright (C) 2007
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: casa-feedback@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: ImagePolProxy.cc 20615 2009-06-09 02:16:01Z Malte.Marquarding $
27 : #include <imageanalysis/ImageAnalysis/ImagePolProxy.h>
28 : #include <casacore/casa/aips.h>
29 : #include <iostream>
30 : #include <sstream>
31 : #include <casacore/casa/IO/ArrayIO.h>
32 : #include <casacore/casa/Arrays/ArrayMath.h>
33 : #include <casacore/casa/Arrays/ArrayUtil.h>
34 : #include <casacore/casa/Arrays/MaskedArray.h>
35 : #include <casacore/casa/Arrays/MaskArrMath.h>
36 : #include <casacore/casa/BasicMath/Random.h>
37 : #include <casacore/casa/BasicSL/String.h>
38 : #include <casacore/casa/Containers/Record.h>
39 : #include <casacore/casa/Exceptions/Error.h>
40 : #include <casacore/casa/Logging/LogFilter.h>
41 : #include <casacore/casa/Logging/LogIO.h>
42 : #include <casacore/casa/Logging/LogOrigin.h>
43 : #include <casacore/coordinates/Coordinates/CoordinateUtil.h>
44 : #include <casacore/coordinates/Coordinates/StokesCoordinate.h>
45 : #include <casacore/coordinates/Coordinates/CoordinateSystem.h>
46 : #include <casacore/coordinates/Coordinates/SpectralCoordinate.h>
47 : #include <casacore/images/Images/ImageExpr.h>
48 : #include <casacore/images/Images/ImageInterface.h>
49 : #include <casacore/images/Regions/ImageRegion.h>
50 : #include <casacore/images/Images/ImageUtilities.h>
51 : #include <casacore/images/Images/PagedImage.h>
52 : #include <casacore/images/Regions/RegionHandler.h>
53 : #include <casacore/images/Images/SubImage.h>
54 : #include <casacore/images/Images/TempImage.h>
55 : #include <casacore/lattices/Lattices/LatticeUtilities.h>
56 : #include <casacore/measures/Measures/Stokes.h>
57 : #include <casacore/tables/Tables/Table.h>
58 : #include <casacore/tables/LogTables/NewFile.h>
59 : #include <casacore/casa/namespace.h>
60 :
61 : #include <memory>
62 :
63 : using namespace casacore;
64 : namespace casa { //# name space casa begins
65 :
66 :
67 0 : ImagePol::ImagePol() : itsImPol(0)
68 : {
69 0 : itsLog= new LogIO();
70 :
71 0 : }
72 :
73 44 : ImagePol::ImagePol(ImageInterface<Float>& im) : itsLog(0),itsImPol(0)
74 : {
75 44 : open(im);
76 :
77 38 : }
78 :
79 76 : ImagePol::~ImagePol(){
80 38 : if(itsLog)
81 38 : delete itsLog;
82 38 : if(itsImPol)
83 38 : delete itsImPol;
84 76 : }
85 :
86 44 : Bool ImagePol::open(ImageInterface<Float>& im){
87 44 : if(!itsLog)
88 44 : itsLog= new LogIO();
89 44 : if(itsImPol)
90 0 : delete itsImPol;
91 44 : itsImPol= new ImagePolarimetry(im);
92 :
93 38 : return true;
94 :
95 : }
96 :
97 0 : Bool ImagePol::open(const String& infile) {
98 0 : Bool rstat(false);
99 :
100 0 : if(!itsLog) itsLog = new LogIO();
101 0 : *itsLog << LogOrigin("ImagePol", "open(const String& infile)");
102 0 : if (infile.size() == 0) {
103 0 : *itsLog << LogIO::WARN << "You must give an infile" << LogIO::POST;
104 0 : return rstat;
105 : }
106 :
107 0 : if(itsImPol) delete itsImPol;
108 0 : itsImPol=0;
109 : //
110 0 : ImageInterface<Float>* imagePointer = 0;
111 0 : ImageUtilities::openImage (imagePointer, infile);
112 :
113 : //
114 : try {
115 0 : itsImPol = new ImagePolarimetry(*imagePointer);
116 0 : rstat = true;
117 0 : } catch (AipsError x) {
118 0 : delete imagePointer;
119 0 : *itsLog << x.getMesg() << LogIO::EXCEPTION;
120 0 : }
121 : //
122 0 : delete imagePointer;
123 0 : imagePointer = 0;
124 0 : return rstat;
125 : }
126 :
127 : Bool
128 0 : ImagePol::imagepoltestimage(const String& outFile, const Vector<Double>& rm,
129 : Bool rmDefault, Double pa0, Double sigma,
130 : Int nx, Int ny, Int nf, Double f0, Double df)
131 : {
132 :
133 0 : Bool rstat(false);
134 :
135 0 : if (itsLog==0) itsLog= new LogIO();
136 0 : *itsLog << LogOrigin("ImagePol", "imagepoltestimage");
137 :
138 0 : if (outFile.size() == 0) {
139 0 : *itsLog << LogIO::WARN << "You must give an outfile" << LogIO::POST;
140 0 : return rstat;
141 : }
142 :
143 0 : if(itsImPol) delete itsImPol;
144 0 : itsImPol = 0;
145 :
146 : // If not given make RM with no ambiguity
147 0 : Vector<Float> rm2(rm.size());
148 0 : if (rmDefault) {
149 0 : Double l1 = QC::c( ).getValue(Unit("m/s")) / f0;
150 0 : Double l2 = QC::c( ).getValue(Unit("m/s")) / (f0+df);
151 0 : rm2.resize(1);
152 0 : rm2(0) = C::pi / 2 / (l1*l1 - l2*l2);
153 : } else {
154 0 : for (uInt i = 0; i < rm.size(); i++) rm2[i] = rm[i];
155 : }
156 :
157 0 : const uInt nRM = rm2.nelements();
158 : //
159 0 : if (nRM == 1) {
160 0 : *itsLog << LogIO::NORMAL << "Using Rotation Measure = " << rm2(0)
161 0 : << " radians/m/m" << endl;
162 : } else {
163 0 : *itsLog << LogIO::NORMAL << "Using Rotation Measures : " << endl;
164 0 : for (uInt i=0; i<nRM; i++) {
165 0 : *itsLog << " " << rm2(i)
166 0 : << " radians/m/m" << endl;
167 : }
168 : }
169 :
170 0 : *itsLog << "Using pa0 = " << pa0 << " degrees" << endl;
171 0 : *itsLog << "Using frequency = " << f0 << " Hz" << endl;
172 0 : *itsLog << "Using bandwidth = " << df << " Hz " << endl;
173 0 : *itsLog << "Using number channels = " << nf << LogIO::POST;
174 :
175 : // Make image
176 0 : IPosition shape(4,nx,ny,4,nf);
177 0 : ImageInterface<Float>* pImOut = 0;
178 0 : makeIQUVImage(pImOut, outFile, sigma, pa0*C::pi/180.0, rm2, shape, f0, df);
179 : try {
180 0 : itsImPol = new ImagePolarimetry(*pImOut);
181 0 : rstat = true;
182 0 : } catch (AipsError x) {
183 0 : delete pImOut;
184 0 : pImOut = 0;
185 0 : *itsLog << x.getMesg() << LogIO::EXCEPTION;
186 0 : }
187 : //
188 0 : delete pImOut;
189 0 : pImOut = 0;
190 :
191 0 : return rstat;
192 0 : }
193 :
194 :
195 0 : Bool ImagePol::depolratio(ImageInterface<Float>*& returnim,
196 : const String& infile, Bool debias,
197 : Double clip, Double sigma,
198 : const String& outfile){
199 0 : Bool rstat(false);
200 0 : *itsLog << LogOrigin("imagepol", __FUNCTION__);
201 0 : if(itsImPol==0){
202 0 : *itsLog <<"No attached image, please use open "
203 0 : << LogIO::EXCEPTION;
204 0 : return rstat;
205 : }
206 0 : auto im1=itsImPol->imageInterface();
207 0 : std::unique_ptr<ImageInterface<Float> > im2;
208 0 : ImageUtilities::openImage(im2, infile);
209 0 : ImageExpr<Float> tmpim=itsImPol->depolarizationRatio(*im1, *im2,
210 : debias,
211 : Float(clip),
212 0 : Float(sigma));
213 :
214 0 : String err;
215 0 : if(!copyImage(returnim, tmpim, outfile, true)){
216 0 : *itsLog <<"Could not convert ratio image "
217 0 : << LogIO::EXCEPTION;
218 : }
219 0 : rstat = true;
220 0 : return rstat;
221 :
222 0 : }
223 :
224 :
225 0 : Bool ImagePol::complexlinpol(const String& outfile) {
226 0 : *itsLog << LogOrigin("ImagePolProxy", __FUNCTION__);
227 0 : if (! itsImPol) {
228 0 : *itsLog << LogIO::SEVERE <<"No attached image, please use open "
229 0 : << LogIO::POST;
230 0 : return false;
231 : }
232 :
233 0 : if (outfile.size() == 0) {
234 0 : *itsLog << LogIO::WARN << "You must give an outfile" << LogIO::POST;
235 0 : return false;
236 : }
237 0 : CoordinateSystem cSysPol;
238 0 : const auto shapePol = itsImPol->singleStokesShape(
239 : cSysPol, Stokes::Plinear
240 0 : );
241 0 : const auto isMasked = itsImPol->isMasked();
242 : auto pOutComplex = _makeImage(
243 : outfile, cSysPol, shapePol, isMasked, false
244 0 : );
245 0 : auto expr = itsImPol->complexLinearPolarization();
246 0 : fiddleStokesCoordinate(*pOutComplex, Stokes::Plinear);
247 :
248 : // Copy to output
249 :
250 0 : pOutComplex->setCoordinateInfo(expr.coordinates());
251 0 : LatticeUtilities::copyDataAndMask(*itsLog, *pOutComplex, expr);
252 0 : auto p = itsImPol->imageInterface();
253 0 : copyMiscellaneous (*pOutComplex, *p);
254 0 : return true;
255 0 : }
256 :
257 : // Summary
258 0 : void ImagePol::summary() const {
259 0 : *itsLog << LogOrigin("imagepol", "summary()");
260 0 : if(itsImPol==0){
261 0 : *itsLog << LogIO::SEVERE <<"No attached image, please use open "
262 0 : << LogIO::POST;
263 0 : return;
264 : }
265 0 : itsImPol->summary(*itsLog);
266 : }
267 :
268 : // sigma
269 0 : Float ImagePol::sigma(Float clip) const {
270 0 : *itsLog << LogOrigin("imagepol", __FUNCTION__);
271 0 : if(itsImPol==0){
272 0 : *itsLog << LogIO::SEVERE <<"No attached image, please use open "
273 0 : << LogIO::POST;
274 0 : return -1.0;
275 : }
276 0 : return itsImPol->sigma(clip);
277 : }
278 :
279 : // Stokes I
280 0 : Bool ImagePol::stokesI(ImageInterface<Float>*& rtnim, const String& outfile){
281 0 : Bool rstat(false);
282 0 : *itsLog << LogOrigin("imagepol", "stokesI(...)");
283 0 : if(itsImPol==0){
284 0 : *itsLog << LogIO::SEVERE <<"No attached image, please use open "
285 0 : << LogIO::POST;
286 0 : return rstat;
287 : }
288 0 : ImageExpr<Float> expr = itsImPol->stokesI();
289 : // Create output image
290 0 : rstat = copyImage (rtnim, expr, outfile, true);
291 0 : return rstat;
292 0 : }
293 :
294 0 : Float ImagePol::sigmaStokesI(Float clip) const {
295 0 : *itsLog << LogOrigin("imagepol", __FUNCTION__);
296 0 : if(itsImPol==0){
297 0 : *itsLog << LogIO::SEVERE <<"No attached image, please use open "
298 0 : << LogIO::POST;
299 0 : return -1.0;
300 : }
301 0 : return itsImPol->sigmaStokesI(clip);
302 : }
303 :
304 : // Stokes Q
305 0 : Bool ImagePol::stokesQ(ImageInterface<Float>*& rtnim, const String& outfile){
306 0 : Bool rstat(false);
307 0 : *itsLog << LogOrigin("imagepol", "stokesQ(...)");
308 0 : if(itsImPol==0){
309 0 : *itsLog << LogIO::SEVERE <<"No attached image, please use open "
310 0 : << LogIO::POST;
311 0 : return rstat;
312 : }
313 0 : ImageExpr<Float> expr = itsImPol->stokesQ();
314 :
315 : // Create output image if needed
316 0 : rstat = copyImage (rtnim, expr, outfile, true);
317 0 : return rstat;
318 0 : }
319 :
320 0 : Float ImagePol::sigmaStokesQ(Float clip) const {
321 0 : *itsLog << LogOrigin("imagepol", "sigmaStokesQ(...)");
322 0 : if(itsImPol==0){
323 0 : *itsLog << LogIO::SEVERE <<"No attached image, please use open "
324 0 : << LogIO::POST;
325 0 : return -1.0;
326 : }
327 0 : return itsImPol->sigmaStokesQ(clip);
328 : }
329 :
330 : // Stokes U
331 0 : Bool ImagePol::stokesU(ImageInterface<Float>*& rtnim, const String& outfile){
332 0 : Bool rstat(false);
333 0 : *itsLog << LogOrigin("imagepol", "stokesU(...)");
334 0 : if(itsImPol==0){
335 0 : *itsLog << LogIO::SEVERE <<"No attached image, please use open "
336 0 : << LogIO::POST;
337 0 : return rstat;
338 : }
339 0 : ImageExpr<Float> expr = itsImPol->stokesU();
340 :
341 : // Create output image if needed
342 0 : rstat = copyImage (rtnim, expr, outfile, true);
343 0 : return rstat;
344 0 : }
345 :
346 0 : Float ImagePol::sigmaStokesU(Float clip) const {
347 0 : *itsLog << LogOrigin("imagepol", "sigmaStokesU(...)");
348 0 : if(itsImPol==0){
349 0 : *itsLog << LogIO::SEVERE <<"No attached image, please use open "
350 0 : << LogIO::POST;
351 0 : return -1.0;
352 : }
353 0 : return itsImPol->sigmaStokesU(clip);
354 : }
355 :
356 : // Stokes V
357 0 : Bool ImagePol::stokesV(ImageInterface<Float>*& rtnim, const String& outfile){
358 0 : Bool rstat(false);
359 0 : *itsLog << LogOrigin("imagepol", "stokesV(...)");
360 0 : if(itsImPol==0){
361 0 : *itsLog << LogIO::SEVERE <<"No attached image, please use open "
362 0 : << LogIO::POST;
363 0 : return rstat;
364 : }
365 0 : ImageExpr<Float> expr = itsImPol->stokesV();
366 :
367 : // Create output image if needed
368 0 : rstat = copyImage (rtnim, expr, outfile, true);
369 0 : return rstat;
370 0 : }
371 :
372 0 : Float ImagePol::sigmaStokesV(Float clip) const {
373 0 : *itsLog << LogOrigin("imagepol", "sigmaStokesV(...)");
374 0 : if(itsImPol==0){
375 0 : *itsLog << LogIO::SEVERE <<"No attached image, please use open "
376 0 : << LogIO::POST;
377 0 : return -1.0;
378 : }
379 0 : return itsImPol->sigmaStokesV(clip);
380 : }
381 :
382 0 : Float ImagePol::sigmaLinPolInt(Float clip, Float sigma) const {
383 0 : *itsLog << LogOrigin("imagepol", __FUNCTION__);
384 0 : if(itsImPol==0){
385 0 : *itsLog << LogIO::SEVERE <<"No attached image, please use open "
386 0 : << LogIO::POST;
387 0 : return -1.0;
388 : }
389 0 : return itsImPol->sigmaLinPolInt(clip, sigma);
390 : }
391 :
392 0 : Float ImagePol::sigmaTotPolInt(Float clip, Float sigma) const {
393 0 : *itsLog << LogOrigin("imagepol", __FUNCTION__);
394 0 : if(itsImPol==0){
395 0 : *itsLog << LogIO::SEVERE <<"No attached image, please use open "
396 0 : << LogIO::POST;
397 0 : return -1.0;
398 : }
399 0 : return itsImPol->sigmaTotPolInt(clip, sigma);
400 : }
401 :
402 : // Complex linear polarization
403 0 : void ImagePol::complexLinearPolarization (const String& outfile) {
404 0 : *itsLog << LogOrigin("imagepol", "complexLinearPolarization(...)");
405 0 : if(itsImPol==0){
406 0 : *itsLog << LogIO::SEVERE <<"No attached image, please use open "
407 0 : << LogIO::POST;
408 0 : return;
409 : }
410 :
411 0 : if (outfile.size() == 0)
412 0 : *itsLog << LogIO::WARN << "You must give an outfile" << LogIO::POST;
413 :
414 : // Make output complex image
415 0 : ImageInterface<Complex>* pOutComplex = 0;
416 0 : CoordinateSystem cSysPol;
417 0 : IPosition shapePol = itsImPol->singleStokesShape(cSysPol, Stokes::Plinear);
418 0 : _makeImage (pOutComplex, outfile, cSysPol, shapePol,
419 0 : itsImPol->isMasked(), false);
420 :
421 : // Make Expr
422 0 : ImageExpr<Complex> expr = itsImPol->complexLinearPolarization();
423 0 : fiddleStokesCoordinate(*pOutComplex, Stokes::Plinear);
424 :
425 : // Copy to output
426 0 : pOutComplex->setCoordinateInfo(expr.coordinates());
427 0 : LatticeUtilities::copyDataAndMask(*itsLog, *pOutComplex, expr);
428 : //
429 0 : auto p = itsImPol->imageInterface();
430 0 : copyMiscellaneous (*pOutComplex, *p);
431 0 : delete pOutComplex;
432 0 : }
433 :
434 : // Complex linear polarization
435 0 : void ImagePol::complexFractionalLinearPolarization (const String& outfile) {
436 0 : *itsLog << LogOrigin("imagepol", "complexFractionalLinearPolarization(...)");
437 0 : if(itsImPol==0){
438 0 : *itsLog << LogIO::SEVERE <<"No attached image, please use open "
439 0 : << LogIO::POST;
440 0 : return;
441 : }
442 :
443 : // Make output complex image
444 0 : ImageInterface<Complex>* pOutComplex = 0;
445 0 : CoordinateSystem cSysPol;
446 0 : IPosition shapePol = itsImPol->singleStokesShape(cSysPol, Stokes::PFlinear);
447 0 : _makeImage (pOutComplex, outfile, cSysPol, shapePol,
448 0 : itsImPol->isMasked(), false);
449 :
450 0 : std::unique_ptr<ImageInterface<Complex> > x(pOutComplex);
451 :
452 : // Make Expr
453 0 : ImageExpr<Complex> expr = itsImPol->complexFractionalLinearPolarization();
454 0 : fiddleStokesCoordinate(*pOutComplex, Stokes::PFlinear);
455 :
456 : // Copy to output
457 0 : pOutComplex->setCoordinateInfo(expr.coordinates());
458 0 : LatticeUtilities::copyDataAndMask(*itsLog, *pOutComplex, expr);
459 0 : copyMiscellaneous (*pOutComplex, *(itsImPol->imageInterface()));
460 0 : }
461 :
462 0 : Bool ImagePol::sigmaLinPolPosAng(ImageInterface<Float>*& rtnim, Float clip,
463 : Float sigma, const String& outfile) {
464 0 : Bool rstat(false);
465 0 : *itsLog << LogOrigin("imagepol", __FUNCTION__);
466 0 : if(itsImPol==0){
467 0 : *itsLog << LogIO::SEVERE <<"No attached image, please use open "
468 0 : << LogIO::POST;
469 0 : return rstat;
470 : }
471 :
472 0 : Bool radians = false;
473 0 : ImageExpr<Float> expr = itsImPol->sigmaLinPolPosAng(radians, clip, sigma);
474 :
475 : // Create output image if needed
476 0 : rstat = copyImage (rtnim, expr, outfile, true);
477 0 : return rstat;
478 0 : }
479 :
480 :
481 : // Fractional linearly polarized intensity
482 0 : Bool ImagePol::fracLinPol(ImageInterface<Float>*& rtnim, Bool debias,
483 : Float clip, Float sigma, const String& outfile) {
484 0 : Bool rstat(false);
485 0 : *itsLog << LogOrigin("imagepol", "fracLinPol(...)");
486 0 : if(itsImPol==0){
487 0 : *itsLog << LogIO::SEVERE <<"No attached image, please use open "
488 0 : << LogIO::POST;
489 0 : return rstat;
490 : }
491 :
492 0 : ImageExpr<Float> expr = itsImPol->fracLinPol(debias, clip, sigma);
493 :
494 : // Create output image
495 0 : rstat = copyImage (rtnim, expr, outfile, true);
496 0 : return rstat;
497 0 : }
498 :
499 :
500 0 : Bool ImagePol::sigmaFracLinPol(ImageInterface<Float>*& rtnim, Float clip,
501 : Float sigma, const String& outfile){
502 0 : Bool rstat(false);
503 0 : *itsLog << LogOrigin("imagepol", __FUNCTION__);
504 0 : if(itsImPol==0){
505 0 : *itsLog << LogIO::SEVERE <<"No attached image, please use open "
506 0 : << LogIO::POST;
507 0 : return rstat;
508 : }
509 :
510 0 : ImageExpr<Float> expr = itsImPol->sigmaFracLinPol(clip, sigma);
511 :
512 : // Create output image if needed
513 0 : rstat = copyImage (rtnim, expr, outfile, true);
514 0 : return rstat;
515 0 : }
516 :
517 :
518 : // Fractional total polarized intensity
519 0 : Bool ImagePol::fracTotPol(ImageInterface<Float>*& rtnim, Bool debias,
520 : Float clip, Float sigma, const String& outfile) {
521 0 : Bool rstat(false);
522 0 : *itsLog << LogOrigin("imagepol", "fracTotPol(...)");
523 0 : if(itsImPol==0){
524 0 : *itsLog << LogIO::SEVERE <<"No attached image, please use open "
525 0 : << LogIO::POST;
526 0 : return rstat;
527 : }
528 :
529 0 : ImageExpr<Float> expr = itsImPol->fracTotPol(debias, clip, sigma);
530 :
531 : // Create output image if needed
532 0 : rstat = copyImage (rtnim, expr, outfile, true);
533 0 : return rstat;
534 0 : }
535 :
536 :
537 0 : Bool ImagePol::sigmaFracTotPol(ImageInterface<Float>*& rtnim, Float clip,
538 : Float sigma, const String& outfile){
539 0 : Bool rstat(false);
540 0 : *itsLog << LogOrigin("imagepol", __FUNCTION__);
541 0 : if(itsImPol==0){
542 0 : *itsLog << LogIO::SEVERE <<"No attached image, please use open "
543 0 : << LogIO::POST;
544 0 : return rstat;
545 : }
546 :
547 0 : ImageExpr<Float> expr = itsImPol->sigmaFracTotPol(clip, sigma);
548 :
549 : // Create output image if needed
550 0 : rstat = copyImage (rtnim, expr, outfile, true);
551 0 : return rstat;
552 0 : }
553 :
554 :
555 : // Depolarization ratio
556 0 : Bool ImagePol::depolarizationRatio (ImageInterface<Float>*& rtnim,
557 : const String& infile,
558 : Bool debias, Float clip,
559 : Float sigma, const String& outfile) {
560 0 : Bool rstat(false);
561 0 : *itsLog << LogOrigin("imagepol", __FUNCTION__);
562 0 : if(itsImPol==0){
563 0 : *itsLog << "No attached image, please use open "
564 0 : << LogIO::EXCEPTION;
565 : }
566 0 : auto imagePointer1 = itsImPol->imageInterface();
567 0 : ImageInterface<Float>* imagePointer2 = 0;
568 0 : ImageUtilities::openImage (imagePointer2, infile);
569 : //
570 : ImageExpr<Float> expr =
571 0 : ImagePolarimetry::depolarizationRatio(*imagePointer1, *imagePointer2,
572 0 : debias, clip, sigma);
573 : //
574 0 : if (imagePointer2) delete imagePointer2;
575 0 : imagePointer2 = 0;
576 :
577 : // Create output image
578 0 : rstat = copyImage (rtnim, expr, outfile, true);
579 :
580 0 : return rstat;
581 0 : }
582 :
583 :
584 0 : Bool ImagePol::sigmaDepolarizationRatio (ImageInterface<Float>*& rtnim,
585 : const String& infile,
586 : Bool debias, Float clip,
587 : Float sigma, const String& outfile) {
588 0 : Bool rstat(false);
589 0 : *itsLog << LogOrigin("imagepol", __FUNCTION__);
590 0 : if(itsImPol==0){
591 0 : *itsLog << LogIO::SEVERE <<"No attached image, please use open "
592 0 : << LogIO::POST;
593 0 : return rstat;
594 : }
595 : //
596 0 : auto imagePointer1 = itsImPol->imageInterface();
597 : //
598 0 : ImageInterface<Float>* imagePointer2 = 0;
599 0 : ImageUtilities::openImage(imagePointer2, infile);
600 : //
601 : ImageExpr<Float> expr =
602 0 : ImagePolarimetry::sigmaDepolarizationRatio(*imagePointer1,
603 : *imagePointer2,
604 0 : debias, clip, sigma);
605 : //
606 0 : if (imagePointer2) delete imagePointer2;
607 0 : imagePointer2 = 0;
608 :
609 : // Create output image
610 0 : rstat = copyImage (rtnim, expr, outfile, true);
611 0 : return rstat;
612 0 : }
613 :
614 :
615 : // Find Rotation Measure from Fourier method
616 0 : void ImagePol::fourierRotationMeasure(const String& outfile,
617 : const String& outfileAmp,
618 : const String& outfilePA,
619 : const String& outfileReal,
620 : const String& outfileImag,
621 : Bool zeroZeroLag) {
622 :
623 0 : *itsLog << LogOrigin("imagepol", __FUNCTION__, WHERE);
624 0 : if(itsImPol==0){
625 0 : *itsLog <<"No attached image, please use open "
626 0 : << LogIO::EXCEPTION;
627 0 : return;
628 : }
629 0 : if (itsImPol->imageInterface()->imageInfo().hasMultipleBeams()) {
630 0 : *itsLog << "Cannot run " << __FUNCTION__
631 : << " on an image with multiple beams"
632 0 : << LogIO::EXCEPTION;
633 : }
634 :
635 :
636 : // Make output complex image
637 0 : ImageInterface<Complex>* pOutComplex = 0;
638 0 : CoordinateSystem cSysPol;
639 0 : IPosition shapePol = itsImPol->singleStokesShape(cSysPol, Stokes::Plinear);
640 0 : _makeImage (pOutComplex, outfile, cSysPol, shapePol,
641 0 : itsImPol->isMasked(), true);
642 :
643 : // Make output amplitude and position angle images
644 0 : ImageInterface<Float>* pOutAmp = 0;
645 0 : ImageInterface<Float>* pOutPA = 0;
646 0 : makeImage (pOutAmp, outfileAmp, cSysPol, shapePol,
647 0 : itsImPol->isMasked(), false);
648 0 : makeImage (pOutPA, outfilePA, cSysPol, shapePol,
649 0 : itsImPol->isMasked(), false);
650 :
651 : // Make output real and imaginary images
652 0 : ImageInterface<Float>* pOutReal = 0;
653 0 : ImageInterface<Float>* pOutImag = 0;
654 0 : makeImage (pOutReal, outfileReal, cSysPol, shapePol,
655 0 : itsImPol->isMasked(), false);
656 0 : makeImage (pOutImag, outfileImag, cSysPol, shapePol,
657 0 : itsImPol->isMasked(), false);
658 :
659 : // The output complex image will have correct Coordinates, mask, and
660 : // miscellaneous things copied to it
661 0 : itsImPol->fourierRotationMeasure(*pOutComplex, zeroZeroLag);
662 :
663 : // Copy to output
664 0 : auto p = itsImPol->imageInterface();
665 0 : if (pOutAmp!=0) {
666 0 : LatticeExprNode node(abs(*pOutComplex));
667 0 : LatticeExpr<Float> le(node);
668 0 : LatticeUtilities::copyDataAndMask(*itsLog, *pOutAmp, le);
669 : //
670 0 : pOutAmp->setCoordinateInfo(pOutComplex->coordinates());
671 0 : copyMiscellaneous (*pOutAmp, *p);
672 0 : pOutAmp->setUnits(p->units());
673 0 : fiddleStokesCoordinate(*pOutAmp, Stokes::Plinear);
674 0 : delete pOutAmp;
675 0 : }
676 0 : if (pOutPA!=0) {
677 0 : LatticeExprNode node(pa(imag(*pOutComplex),real(*pOutComplex))); // degrees
678 0 : LatticeExpr<Float> le(node);
679 0 : LatticeUtilities::copyDataAndMask(*itsLog, *pOutPA, le);
680 : //
681 0 : pOutPA->setCoordinateInfo(pOutComplex->coordinates());
682 0 : copyMiscellaneous (*pOutPA, *p);
683 0 : pOutPA->setUnits("deg");
684 0 : fiddleStokesCoordinate(*pOutPA, Stokes::Pangle);
685 0 : delete pOutPA;
686 0 : }
687 0 : if (pOutReal!=0) {
688 0 : LatticeExprNode node(real(*pOutComplex));
689 0 : LatticeExpr<Float> le(node);
690 0 : LatticeUtilities::copyDataAndMask(*itsLog, *pOutReal, le);
691 0 : pOutReal->setCoordinateInfo(pOutComplex->coordinates());
692 0 : copyMiscellaneous (*pOutReal, *p);
693 0 : pOutReal->setUnits(p->units());
694 0 : fiddleStokesCoordinate(*pOutReal, Stokes::Plinear); // Not strictly correct
695 0 : delete pOutReal;
696 0 : }
697 0 : if (pOutImag!=0) {
698 0 : LatticeExprNode node(imag(*pOutComplex));
699 0 : LatticeExpr<Float> le(node);
700 0 : LatticeUtilities::copyDataAndMask(*itsLog, *pOutImag, le);
701 0 : pOutImag->setCoordinateInfo(pOutComplex->coordinates());
702 0 : copyMiscellaneous (*pOutImag, *p);
703 0 : pOutImag->setUnits(p->units());
704 0 : fiddleStokesCoordinate(*pOutImag, Stokes::Plinear); // Not strictly correct
705 0 : delete pOutImag;
706 0 : }
707 0 : }
708 :
709 16 : void ImagePol::rotationMeasure(
710 : const String& outRM, const String& outRMErr,
711 : const String& outPA0, const String& outPA0Err,
712 : const String& outNTurns, const String& outChiSq,
713 : Int axis2, Float sigmaQU, Float rmFg,
714 : Float rmMax, Float maxPaErr
715 : ) {
716 : // Find Rotation Measure from traditional method
717 16 : *itsLog << LogOrigin("ImagePol", __func__);
718 16 : if(itsImPol==0) {
719 0 : *itsLog << LogIO::SEVERE <<"No attached image, please use open "
720 0 : << LogIO::POST;
721 0 : return;
722 : }
723 : // Make output images. Give them all a mask as we don't know if output
724 : // will be masked or not.
725 16 : CoordinateSystem cSysRM;
726 : Int fAxis, sAxis;
727 16 : Int axis = axis2;
728 16 : IPosition shapeRM = itsImPol->rotationMeasureShape(
729 16 : cSysRM, fAxis, sAxis, *itsLog, axis
730 16 : );
731 16 : ImageInterface<Float>* pRMOut = 0;
732 16 : ImageInterface<Float>* pRMOutErr = 0;
733 16 : makeImage (pRMOut, outRM, cSysRM, shapeRM, true, false);
734 : // manage naked pointers so exception throwing doesn't leave open images
735 16 : std::unique_ptr<ImageInterface<Float> > managed5(pRMOut);
736 16 : makeImage (pRMOutErr, outRMErr, cSysRM, shapeRM, true, false);
737 16 : std::unique_ptr<ImageInterface<Float> > managed6(pRMOutErr);
738 16 : CoordinateSystem cSysPA;
739 16 : IPosition shapePA = itsImPol->positionAngleShape(
740 16 : cSysPA, fAxis, sAxis, *itsLog, axis
741 16 : );
742 16 : ImageInterface<Float>* pPA0Out = 0;
743 16 : ImageInterface<Float>* pPA0OutErr = 0;
744 16 : makeImage (pPA0Out, outPA0, cSysPA, shapePA, true, false);
745 16 : std::unique_ptr<ImageInterface<Float> > managed1(pPA0Out);
746 16 : makeImage (pPA0OutErr, outPA0Err, cSysPA, shapePA, true, false);
747 16 : std::unique_ptr<ImageInterface<Float> > managed2(pPA0OutErr);
748 16 : ImageInterface<Float>* pNTurnsOut = 0;
749 16 : makeImage (pNTurnsOut, outNTurns, cSysRM, shapeRM, true, false);
750 16 : std::unique_ptr<ImageInterface<Float> > managed3(pNTurnsOut);
751 16 : ImageInterface<Float>* pChiSqOut = 0;
752 16 : makeImage (pChiSqOut, outChiSq, cSysRM, shapeRM, true, false);
753 16 : std::unique_ptr<ImageInterface<Float> > managed4(pChiSqOut);
754 16 : itsImPol->rotationMeasure(
755 : pRMOut, pRMOutErr, pPA0Out, pPA0OutErr,
756 : pNTurnsOut, pChiSqOut,
757 : axis, rmMax, maxPaErr,
758 : sigmaQU, rmFg, true
759 : );
760 15 : auto p = itsImPol->imageInterface();
761 15 : if (pRMOut) {
762 15 : copyMiscellaneous (*pRMOut, *p);
763 : }
764 15 : if (pRMOutErr) {
765 0 : copyMiscellaneous (*pRMOutErr, *p);
766 : }
767 15 : if (pPA0Out) {
768 1 : copyMiscellaneous (*pPA0Out, *p);
769 : }
770 15 : if (pPA0OutErr) {
771 0 : copyMiscellaneous (*pPA0OutErr, *p);
772 : }
773 15 : if (pNTurnsOut) {
774 0 : copyMiscellaneous (*pNTurnsOut, *p);
775 : }
776 15 : if (pChiSqOut) {
777 0 : copyMiscellaneous (*pChiSqOut, *p);
778 : }
779 25 : }
780 :
781 :
782 : // Make a complex image
783 0 : void ImagePol::makeComplex (const String& outfile, const String& real,
784 : const String& imag, const String& amp,
785 : const String& phase) {
786 :
787 0 : *itsLog << LogOrigin("imagepol", __FUNCTION__);
788 0 : if(itsImPol==0){
789 0 : *itsLog << LogIO::SEVERE <<"No attached image, please use open "
790 0 : << LogIO::POST;
791 0 : return;
792 : }
793 :
794 : // Checks
795 0 : if (outfile.empty()) {
796 0 : *itsLog << "You must give the output complex image file name"
797 0 : << LogIO::EXCEPTION;
798 : }
799 :
800 0 : Bool doRI = !real.empty() && !imag.empty();
801 0 : Bool doAP = !amp.empty() && !phase.empty();
802 0 : if (doRI && doAP) {
803 0 : *itsLog << "You must give either real and imaginary, or amplitude and phase; Not all four of them"
804 0 : << LogIO::EXCEPTION;
805 : }
806 :
807 : // Make output complex image
808 0 : ImageInterface<Complex>* pOutComplex = 0;
809 0 : CoordinateSystem cSysPol;
810 0 : IPosition shapePol = itsImPol->singleStokesShape(cSysPol, Stokes::I);
811 0 : _makeImage (pOutComplex, outfile, cSysPol, shapePol,
812 0 : itsImPol->isMasked(), false);
813 :
814 : // Make Expression. Only limited Stokes types that make sense are allowed.
815 0 : ImageExpr<Complex>* pExpr = 0;
816 0 : if (doRI) {
817 0 : PagedImage<Float> rr(real);
818 0 : Stokes::StokesTypes tR = stokesType(rr.coordinates());
819 : //
820 0 : PagedImage<Float> ii(imag);
821 0 : Stokes::StokesTypes tI = stokesType(ii.coordinates());
822 : //
823 0 : if (tR!=Stokes::Q || tI!=Stokes::U) {
824 0 : *itsLog << "The real and imaginary components must be Q and U, respectively"
825 0 : << LogIO::EXCEPTION;
826 : }
827 0 : Stokes::StokesTypes typeOut = Stokes::Plinear;
828 : //
829 0 : LatticeExprNode node(formComplex(rr,ii));
830 0 : LatticeExpr<Complex> le(node);
831 0 : pExpr = new ImageExpr<Complex>(le, String("ComplexLinearPolarization"));
832 0 : fiddleStokesCoordinate(*pExpr, typeOut);
833 0 : } else {
834 0 : PagedImage<Float> aa(amp);
835 0 : Stokes::StokesTypes tA = stokesType(aa.coordinates());
836 : //
837 0 : PagedImage<Float> pp(phase);
838 0 : Stokes::StokesTypes tP = stokesType(pp.coordinates());
839 : //
840 0 : if (tP!=Stokes::Pangle) {
841 0 : *itsLog << "The phase must be of Stokes type position angle (Pangle)"
842 0 : << LogIO::EXCEPTION;
843 : }
844 0 : Float fac = 1.0;
845 0 : String units = pp.units().getName();
846 0 : if (units.contains(String("deg"))) {
847 0 : fac = C::pi / 180.0;
848 0 : } else if (units.contains(String("rad"))) {
849 : } else {
850 0 : *itsLog << LogIO::WARN
851 : << "Units for phase are neither radians nor degrees. radians assumed"
852 0 : << LogIO::POST;
853 : }
854 : //
855 0 : Stokes::StokesTypes typeOut = Stokes::Undefined;
856 0 : String exprName("");
857 0 : if (tA==Stokes::Ptotal) {
858 0 : typeOut = Stokes::Ptotal;
859 0 : exprName = String("ComplexTotalPolarization");
860 0 : } else if (tA==Stokes::Plinear) {
861 0 : typeOut = Stokes::Plinear;
862 0 : exprName = String("ComplexLinearPolarization");
863 0 : } else if (tA==Stokes::PFtotal) {
864 0 : typeOut = Stokes::PFtotal;
865 0 : exprName = String("ComplexFractionalTotalPolarization");
866 0 : } else if (tA==Stokes::PFlinear) {
867 0 : typeOut = Stokes::PFlinear;
868 0 : exprName = String("ComplexFractionalLinearPolarization");
869 : } else {
870 0 : *itsLog << "Cannot form Complex image for this amplitude image"
871 0 : << endl;
872 0 : *itsLog << "Expect linear, total, or fractional polarization"
873 0 : << LogIO::EXCEPTION;
874 : }
875 : //
876 0 : LatticeExprNode node0(2.0*fac*pp);
877 0 : LatticeExprNode node(formComplex(aa*cos(node0),aa*sin(node0)));
878 0 : LatticeExpr<Complex> le(node);
879 0 : pExpr = new ImageExpr<Complex>(le, exprName);
880 0 : fiddleStokesCoordinate(*pExpr, typeOut);
881 0 : }
882 :
883 : // Copy to output
884 0 : pOutComplex->setCoordinateInfo(pExpr->coordinates());
885 0 : LatticeUtilities::copyDataAndMask(*itsLog, *pOutComplex, *pExpr);
886 : //
887 0 : auto p = itsImPol->imageInterface();
888 0 : copyMiscellaneous (*pOutComplex, *p);
889 : //
890 0 : delete pExpr; pExpr = 0;
891 0 : delete pOutComplex; pOutComplex = 0;
892 0 : }
893 :
894 :
895 : // Private functions
896 :
897 0 : Bool ImagePol::copyImage(ImageInterface<Float>*& out,
898 : const ImageInterface<Float>& in,
899 : const String& outfile, Bool overwrite){
900 :
901 : // if no outfile, just make out image from the input image
902 0 : if (outfile.empty()){
903 0 : out = in.cloneII();
904 0 : return true;
905 : }
906 :
907 : // The user wants to write the image out; verify file
908 0 : if (!overwrite) {
909 0 : NewFile validfile;
910 0 : String errmsg;
911 0 : if (!validfile.valueOK(outfile, errmsg)) {
912 0 : *itsLog << errmsg << LogIO::EXCEPTION;
913 : }
914 0 : }
915 :
916 : // Create output image
917 0 : out= new PagedImage<Float>(in.shape(), in.coordinates(), outfile);
918 0 : if (out == 0) {
919 0 : *itsLog << "Failed to create PagedImage" << LogIO::EXCEPTION;
920 : } else {
921 0 : *itsLog << LogIO::NORMAL << "Creating image '" << outfile
922 0 : << "' of shape " << out->shape() << LogIO::POST;
923 : }
924 :
925 : // Make mask
926 0 : if (in.isMasked()) makeMask(*out, false);
927 :
928 : // Copy stuff
929 0 : LatticeUtilities::copyDataAndMask (*itsLog, *out, in);
930 0 : ImageUtilities::copyMiscellaneous (*out, in);
931 :
932 0 : return true;
933 :
934 : }
935 :
936 17 : Bool ImagePol::makeMask(ImageInterface<Float>& out, Bool init){
937 :
938 17 : if (out.canDefineRegion()) {
939 :
940 : // Generate mask name if not given
941 17 : String maskName = out.makeUniqueRegionName(String("mask"), 0);
942 :
943 : // Make the mask if it does not exist
944 17 : if (!out.hasRegion (maskName, RegionHandler::Masks)) {
945 17 : out.makeMask (maskName, true, true, init, true);
946 : /* if (init) {
947 : *itsLog << LogIO::NORMAL << "Created and initialized mask `" << maskName << "'" << LogIO::POST;
948 : } else {
949 : *itsLog << LogIO::NORMAL << "Created mask `" << maskName << "'" << LogIO::POST;
950 : }*/
951 : }
952 17 : return true;
953 17 : } else {
954 0 : *itsLog << LogIO::WARN
955 0 : << "Cannot make requested mask for this type of image" << endl;
956 0 : return false;
957 : }
958 :
959 : }
960 0 : Bool ImagePol::makeMask(ImageInterface<Complex>& out, Bool init){
961 :
962 0 : if (out.canDefineRegion()) {
963 :
964 : // Generate mask name if not given
965 0 : String maskName = out.makeUniqueRegionName(String("mask"), 0);
966 :
967 : // Make the mask if it does not exist
968 0 : if (!out.hasRegion (maskName, RegionHandler::Masks)) {
969 0 : out.makeMask (maskName, true, true, init, true);
970 : /* if (init) {
971 : *itsLog << LogIO::NORMAL << "Created and initialized mask `" << maskName << "'" << LogIO::POST;
972 : } else {
973 : *itsLog << LogIO::NORMAL << "Created mask `" << maskName << "'" << LogIO::POST;
974 : }*/
975 : }
976 0 : return true;
977 0 : } else {
978 0 : *itsLog << LogIO::WARN
979 0 : << "Cannot make requested mask for this type of image" << endl;
980 0 : return false;
981 : }
982 :
983 : }
984 :
985 :
986 96 : Bool ImagePol::makeImage(ImageInterface<Float>*& out,
987 : const String& outfile, const CoordinateSystem& cSys,
988 : const IPosition& shape, Bool isMasked,
989 : Bool tempAllowed) {
990 : // Verify outfile
991 96 : if (outfile.empty()) {
992 79 : if (!tempAllowed) return false;
993 : }
994 : // else {
995 : // NewFile validfile;
996 : // String errmsg;
997 : // if (!validfile.valueOK(outfile, errmsg)) {
998 : // *itsLog << errmsg << LogIO::EXCEPTION;
999 : // }
1000 : //}
1001 :
1002 :
1003 17 : uInt ndim = shape.nelements();
1004 17 : if (ndim != cSys.nPixelAxes()) {
1005 0 : *itsLog << "Supplied CoordinateSystem and image shape are inconsistent"
1006 0 : << LogIO::EXCEPTION;
1007 : }
1008 17 : if (outfile.empty()) {
1009 0 : out = new TempImage<Float>(shape, cSys);
1010 0 : if (out == 0) {
1011 0 : *itsLog << "Failed to create TempImage" << LogIO::EXCEPTION;
1012 : }
1013 0 : *itsLog << LogIO::NORMAL << "Creating (temp)image of shape "
1014 0 : << out->shape() << LogIO::POST;
1015 : } else {
1016 17 : out = new PagedImage<Float>(shape, cSys, outfile);
1017 17 : if (out == 0) {
1018 0 : *itsLog << "Failed to create PagedImage" << LogIO::EXCEPTION;
1019 : }
1020 17 : *itsLog << LogIO::NORMAL << "Creating image '" << outfile
1021 17 : << "' of shape " << out->shape() << LogIO::POST;
1022 : }
1023 : //
1024 17 : if (isMasked) {
1025 17 : makeMask(*out, true);
1026 : }
1027 :
1028 17 : return true;
1029 : }
1030 :
1031 0 : Bool ImagePol::_makeImage(
1032 : ImageInterface<Complex>*& out, const String& outfile,
1033 : const CoordinateSystem& cSys, const IPosition& shape, Bool isMasked,
1034 : Bool tempAllowed
1035 : ) {
1036 0 : const auto ndim = shape.nelements();
1037 0 : if (ndim != cSys.nPixelAxes()) {
1038 0 : *itsLog << "Supplied CoordinateSystem and image shape are inconsistent"
1039 0 : << LogIO::EXCEPTION;
1040 : }
1041 0 : if (outfile.empty()) {
1042 0 : if (tempAllowed) {
1043 0 : out = new TempImage<Complex>(shape, cSys);
1044 0 : if (! out) {
1045 0 : *itsLog << "Failed to create TempImage" << LogIO::EXCEPTION;
1046 : }
1047 0 : *itsLog << LogIO::NORMAL << "Creating (temp)image of shape "
1048 0 : << out->shape() << LogIO::POST;
1049 : }
1050 : else {
1051 0 : return false;
1052 : }
1053 : }
1054 : else {
1055 0 : out = new PagedImage<Complex>(shape, cSys, outfile);
1056 0 : if (! out) {
1057 0 : *itsLog << "Failed to create PagedImage" << LogIO::EXCEPTION;
1058 : }
1059 0 : *itsLog << LogIO::NORMAL << "Creating image '" << outfile
1060 0 : << "' of shape " << out->shape() << LogIO::POST;
1061 : }
1062 0 : if (isMasked) {
1063 0 : makeMask(*out, true);
1064 : }
1065 0 : return true;
1066 : }
1067 :
1068 0 : SPIIC ImagePol::_makeImage(
1069 : const String& outfile, const CoordinateSystem& cSys,
1070 : const IPosition& shape, bool isMasked,
1071 : bool tempAllowed
1072 : ) {
1073 0 : ImageInterface<Complex>* out = NULL;
1074 0 : const auto res = _makeImage(
1075 : out, outfile, cSys, shape, isMasked, tempAllowed
1076 : );
1077 0 : if (res) {
1078 0 : return shared_ptr<ImageInterface<Complex>>(out);
1079 : }
1080 : else {
1081 0 : if (out) {
1082 0 : delete out;
1083 : }
1084 0 : return shared_ptr<ImageInterface<Complex>>();
1085 : }
1086 : }
1087 :
1088 0 : void ImagePol::copyMiscellaneous (ImageInterface<Complex>& out,
1089 : const ImageInterface<Float>& in)
1090 : {
1091 0 : out.setMiscInfo(in.miscInfo());
1092 0 : out.appendLog(in.logger());
1093 0 : }
1094 :
1095 16 : void ImagePol::copyMiscellaneous (ImageInterface<Float>& out,
1096 : const ImageInterface<Float>& in)
1097 : {
1098 16 : out.setMiscInfo(in.miscInfo());
1099 16 : out.appendLog(in.logger());
1100 16 : }
1101 :
1102 :
1103 0 : void ImagePol::fiddleStokesCoordinate(ImageInterface<Float>& ie,
1104 : Stokes::StokesTypes type)
1105 : {
1106 0 : CoordinateSystem cSys = ie.coordinates();
1107 : //
1108 0 : Int afterCoord = -1;
1109 0 : Int iStokes = cSys.findCoordinate(Coordinate::STOKES, afterCoord);
1110 : //
1111 0 : Vector<Int> which(1);
1112 0 : which(0) = Int(type);
1113 0 : StokesCoordinate stokes(which);
1114 0 : cSys.replaceCoordinate(stokes, iStokes);
1115 0 : ie.setCoordinateInfo(cSys);
1116 0 : }
1117 :
1118 0 : void ImagePol::fiddleStokesCoordinate(ImageInterface<Complex>& ie,
1119 : Stokes::StokesTypes type)
1120 : {
1121 0 : CoordinateSystem cSys = ie.coordinates();
1122 : //
1123 0 : Int afterCoord = -1;
1124 0 : Int iStokes = cSys.findCoordinate(Coordinate::STOKES, afterCoord);
1125 : //
1126 0 : Vector<Int> which(1);
1127 0 : which(0) = Int(type);
1128 0 : StokesCoordinate stokes(which);
1129 0 : cSys.replaceCoordinate(stokes, iStokes);
1130 0 : ie.setCoordinateInfo(cSys);
1131 0 : }
1132 :
1133 :
1134 0 : Bool ImagePol::makeIQUVImage (ImageInterface<Float>*& pImOut,
1135 : const String& outfile,
1136 : Double sigma, Double pa0,
1137 : const Vector<Float>& rm,
1138 : const IPosition& shape,
1139 : Double f0, Double dF)
1140 : //
1141 : // Must be 4D
1142 : //
1143 : {
1144 0 : AlwaysAssert(shape.nelements()==4,AipsError);
1145 0 : AlwaysAssert(shape(2)==4,AipsError);
1146 : //
1147 0 : CoordinateSystem cSys;
1148 0 : CoordinateUtil::addDirAxes(cSys);
1149 : //
1150 0 : Vector<Int> whichStokes(4);
1151 0 : whichStokes(0) = Stokes::I;
1152 0 : whichStokes(1) = Stokes::Q;
1153 0 : whichStokes(2) = Stokes::U;
1154 0 : whichStokes(3) = Stokes::V;
1155 0 : StokesCoordinate stokesCoord(whichStokes);
1156 0 : cSys.addCoordinate(stokesCoord);
1157 : //
1158 0 : const Int nchan = shape(3);
1159 0 : Double df = dF / nchan;
1160 0 : Double refpix = 0.0;
1161 0 : SpectralCoordinate spectCoord(MFrequency::TOPO, f0, df, refpix, f0);
1162 0 : cSys.addCoordinate(spectCoord);
1163 :
1164 : // Centre reference pixel
1165 0 : centreRefPix (cSys, shape);
1166 :
1167 : // Make image
1168 0 : makeImage (pImOut, outfile, cSys, shape, false, true);
1169 : //
1170 0 : uInt stokesAxis = 2;
1171 0 : uInt spectralAxis = 3;
1172 :
1173 : // Fill image with I, Q, U and V.
1174 0 : fillIQUV (*pImOut, stokesAxis, spectralAxis, rm, pa0);
1175 :
1176 : // Add noise
1177 0 : Array<Float> slice = pImOut->get();
1178 0 : Float maxVal = max(slice);
1179 0 : Float t = sigma * maxVal;
1180 0 : *itsLog << LogIO::NORMAL << "Using sigma = "
1181 0 : << t << LogIO::POST;
1182 0 : MLCG gen;
1183 0 : Normal noiseGen(&gen, 0.0, t*t);
1184 0 : addNoise(slice, noiseGen);
1185 0 : pImOut->put(slice);
1186 0 : return true;
1187 0 : }
1188 :
1189 0 : Bool ImagePol::fillIQUV (ImageInterface<Float>& im, uInt stokesAxis,
1190 : uInt spectralAxis, const Vector<Float>& rm,
1191 : Float pa0)
1192 : //
1193 : // Image must be 4D
1194 : //
1195 : {
1196 :
1197 : // Find spectral coordinate
1198 0 : const CoordinateSystem& cSys = im.coordinates();
1199 : Int spectralCoord, iDum;
1200 0 : cSys.findPixelAxis(spectralCoord, iDum, spectralAxis);
1201 0 : const SpectralCoordinate& sC = cSys.spectralCoordinate(spectralCoord);
1202 : //
1203 0 : IPosition shape = im.shape();
1204 0 : Double c = QC::c( ).getValue(Unit("m/s"));
1205 : Double lambdasq;
1206 0 : MFrequency freq;
1207 0 : IPosition blc(4,0);
1208 0 : IPosition trc(shape-1);
1209 : //
1210 0 : Double ii = 2.0; // arbitrary
1211 0 : Double vv = 0.05 * ii; // arbitrary
1212 0 : const Int n = shape(3);
1213 0 : const uInt nRM = rm.nelements();
1214 0 : for (Int i=0; i<n; i++) {
1215 0 : if (!sC.toWorld(freq, Double(i))) {
1216 0 : *itsLog << sC.errorMessage() << LogIO::EXCEPTION;
1217 : }
1218 0 : Double fac = c / freq.get(Unit("Hz")).getValue();
1219 0 : lambdasq = fac*fac;
1220 : //
1221 0 : Double chi = rm(0)*lambdasq + pa0;
1222 0 : Double q = cos(2*chi);
1223 0 : Double u = sin(2*chi);
1224 : //
1225 0 : if (nRM > 1) {
1226 0 : for (uInt j=1; j<nRM; j++) {
1227 0 : chi = rm(j)*lambdasq + pa0;
1228 0 : q += cos(2*chi);
1229 0 : u += sin(2*chi);
1230 : }
1231 : }
1232 0 : q = q / Double(nRM);
1233 0 : u = u / Double(nRM);
1234 : //
1235 0 : blc(spectralAxis) = i; // channel
1236 0 : trc(spectralAxis) = i;
1237 : {
1238 0 : blc(stokesAxis) = 1; // Q
1239 0 : trc(stokesAxis) = 1;
1240 0 : Slicer sl(blc, trc, Slicer::endIsLast);
1241 0 : SubImage<Float> subImage(im, sl, true);
1242 0 : subImage.set(q);
1243 0 : }
1244 : {
1245 0 : blc(stokesAxis) = 2; // U
1246 0 : trc(stokesAxis) = 2;
1247 0 : Slicer sl(blc, trc, Slicer::endIsLast);
1248 0 : SubImage<Float> subImage(im, sl, true);
1249 0 : subImage.set(u);
1250 0 : }
1251 : }
1252 : //
1253 0 : blc(spectralAxis) = 0;
1254 0 : trc(spectralAxis) = n-1;
1255 : {
1256 0 : blc(stokesAxis) = 0; // I
1257 0 : trc(stokesAxis) = 0;
1258 0 : Slicer sl(blc, trc, Slicer::endIsLast);
1259 0 : SubImage<Float> subImage(im, sl, true);
1260 0 : subImage.set(ii);
1261 0 : }
1262 : {
1263 0 : blc(stokesAxis) = 3; // V
1264 0 : trc(stokesAxis) = 3;
1265 0 : Slicer sl(blc, trc, Slicer::endIsLast);
1266 0 : SubImage<Float> subImage(im, sl, true);
1267 0 : subImage.set(vv);
1268 0 : }
1269 0 : return true;
1270 :
1271 0 : }
1272 :
1273 0 : void ImagePol::addNoise (Array<Float>& slice, Normal& noiseGen)
1274 : {
1275 : Bool deleteIt;
1276 0 : Float* p = slice.getStorage(deleteIt);
1277 0 : for (uInt i=0; i<slice.nelements(); i++) {
1278 0 : p[i] += noiseGen();
1279 : }
1280 0 : slice.putStorage(p, deleteIt);
1281 0 : }
1282 :
1283 0 : void ImagePol::centreRefPix (CoordinateSystem& cSys, const IPosition& shape)
1284 : {
1285 0 : Int after = -1;
1286 0 : Int iS = cSys.findCoordinate(Coordinate::STOKES, after);
1287 0 : Int sP = -1;
1288 0 : if (iS>=0) {
1289 0 : Vector<Int> pixelAxes = cSys.pixelAxes(iS);
1290 0 : sP = pixelAxes(0);
1291 0 : }
1292 0 : Vector<Double> refPix = cSys.referencePixel();
1293 0 : for (Int i=0; i<Int(refPix.nelements()); i++) {
1294 0 : if (i!=sP) refPix(i) = Double(shape(i) / 2);
1295 : }
1296 0 : cSys.setReferencePixel(refPix);
1297 0 : }
1298 :
1299 :
1300 0 : Stokes::StokesTypes ImagePol::stokesType(const CoordinateSystem& cSys)
1301 : {
1302 0 : Stokes::StokesTypes type = Stokes::Undefined;
1303 0 : Int afterCoord = -1;
1304 0 : Int iStokes = cSys.findCoordinate(Coordinate::STOKES, afterCoord);
1305 0 : if (iStokes >=0) {
1306 0 : Vector<Int> which = cSys.stokesCoordinate(iStokes).stokes();
1307 0 : if (which.nelements()>1) {
1308 0 : *itsLog << "Stokes axis must be of length unity" << LogIO::EXCEPTION;
1309 : } else {
1310 0 : type = Stokes::type(which(0));
1311 : }
1312 0 : } else {
1313 0 : *itsLog << "No StokesCoordinate" << LogIO::EXCEPTION;
1314 : }
1315 0 : return type;
1316 : }
1317 :
1318 :
1319 : } // end of casa namespace
|
