Line data Source code
1 : //# StokesImageUtil.cc: Stokes Image Utilities
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 adressed 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 : //#
27 : //# $Id$
28 :
29 : #include <casacore/casa/aips.h>
30 :
31 : #include <casacore/casa/Arrays/Matrix.h>
32 : #include <casacore/casa/Arrays/MatrixMath.h>
33 : #include <casacore/casa/Arrays/ArrayMath.h>
34 : #include <casacore/casa/Arrays/Cube.h>
35 : #include <casacore/casa/BasicMath/Math.h>
36 : #include <casacore/casa/BasicSL/Complex.h>
37 : #include <casacore/scimath/Mathematics/FFTServer.h>
38 : #include <casacore/casa/BasicSL/Constants.h>
39 : #include <casacore/casa/Utilities/Assert.h>
40 : #include <casacore/lattices/Lattices/Lattice.h>
41 : #include <casacore/lattices/LEL/LatticeExpr.h>
42 : #include <casacore/lattices/Lattices/LatticeIterator.h>
43 : #include <casacore/lattices/Lattices/LatticeStepper.h>
44 : #include <casacore/lattices/LatticeMath/LatticeConvolver.h>
45 : #include <casacore/scimath/Fitting/NonLinearFitLM.h>
46 : #include <casacore/scimath/Functionals/Gaussian2D.h>
47 : #include <casacore/scimath/Mathematics/Interpolate2D.h>
48 : #include <casacore/casa/IO/ArrayIO.h>
49 :
50 : #include <casacore/ms/MeasurementSets/MSColumns.h>
51 : #include <casacore/ms/MeasurementSets/MSDopplerUtil.h>
52 :
53 : #include <msvis/MSVis/StokesVector.h>
54 : #include <synthesis/TransformMachines/StokesImageUtil.h>
55 : #include <casacore/images/Images/PagedImage.h>
56 : #include <casacore/images/Images/SubImage.h>
57 : #include <casacore/images/Images/ImageBeamSet.h>
58 : #include <casacore/casa/OS/File.h>
59 :
60 : #include <casacore/casa/Quanta/UnitMap.h>
61 : #include <casacore/casa/Quanta/UnitVal.h>
62 : #include <casacore/measures/Measures/Stokes.h>
63 : #include <casacore/coordinates/Coordinates/CoordinateSystem.h>
64 : #include <casacore/coordinates/Coordinates/DirectionCoordinate.h>
65 : #include <casacore/coordinates/Coordinates/SpectralCoordinate.h>
66 : #include <casacore/coordinates/Coordinates/StokesCoordinate.h>
67 : #include <casacore/coordinates/Coordinates/Projection.h>
68 : #include <casacore/casa/Logging/LogSink.h>
69 : #include <casacore/casa/Logging/LogMessage.h>
70 :
71 : #include <synthesis/TransformMachines2/Utils.h>
72 :
73 : #include <iostream>
74 : #include <ctime>
75 :
76 : using namespace casacore;
77 : namespace casa {
78 :
79 :
80 :
81 : // <summary>
82 : // </summary>
83 :
84 : // <reviewed reviewer="" date="" tests="tMEGI" demos="">
85 :
86 : // <prerequisite>
87 : // </prerequisite>
88 : //
89 : // <etymology>
90 : // </etymology>
91 : //
92 : // <synopsis>
93 : // </synopsis>
94 : //
95 : // <example>
96 : // <srcblock>
97 : // </srcblock>
98 : // </example>
99 : //
100 : // <motivation>
101 : // </motivation>
102 : //
103 : // <todo asof="">
104 : // </todo>
105 :
106 : // In-place convolve. image is 4 dimensional. RA and Dec first.
107 3130 : void StokesImageUtil::Convolve(ImageInterface<Float>& image,
108 : ImageInterface<Float>& psf) {
109 :
110 :
111 3130 : Double availablemem=Double(HostInfo::memoryFree())*1024.0;
112 : //Do an in memory if x-y complex can fit in remaining memory
113 3130 : if(availablemem < Double(image.shape()(0)*image.shape()(1))*6.0*8.0 ){
114 0 : LatticeConvolver<Float> lattCon(psf);
115 0 : lattCon.linear(image);
116 0 : }
117 : else{
118 :
119 3130 : Vector<Int> map;
120 3130 : AlwaysAssert(StokesMap(map, image.coordinates()), AipsError);
121 :
122 3130 : Vector<Int> psfmap;
123 3130 : AlwaysAssert(StokesMap(psfmap, psf.coordinates()), AipsError);
124 :
125 3130 : Int nx = image.shape()(map(0));
126 3130 : Int ny = image.shape()(map(1));
127 :
128 : // Make FFT machine
129 3130 : FFTServer<Float,Complex> fft(IPosition(2, nx, ny));
130 :
131 : // Get the image into a Matrix for Fourier transformation
132 3130 : Array<Complex> xfr;
133 3130 : Matrix<Complex> cft;
134 :
135 6260 : LatticeStepper ls(image.shape(), IPosition(4, nx, ny, 1, 1),
136 6260 : IPosition(4, map(0), map(1), map(2), map(3)));
137 3130 : LatticeIterator<Float> li(image, ls);
138 :
139 : // Get the PSF into a Matrix for Fourier transformation
140 6260 : LatticeStepper psfls(psf.shape(), IPosition(4, nx, ny, 1, 1),
141 6260 : IPosition(4, map(0), map(1), map(2), map(3)));
142 3130 : RO_LatticeIterator<Float> psfli(psf, psfls);
143 3130 : psfli.reset();
144 : //fft.fft(xfr, psfli.matrixCursor());
145 3130 : fft.fft0(xfr, psfli.matrixCursor());
146 :
147 : // Loop over all planes
148 6260 : for (li.reset();!li.atEnd();li++) {
149 : //fft.fft(cft, li.matrixCursor());
150 3130 : fft.fft0(cft,li.matrixCursor());
151 3130 : cft *= xfr;
152 : // fft.fft(li.rwMatrixCursor(), cft);
153 3130 : fft.fft0(li.rwMatrixCursor(),cft,false);
154 3130 : fft.flip(li.rwMatrixCursor(), false, false);
155 : }
156 3130 : }
157 3130 : };
158 :
159 0 : void StokesImageUtil::Convolve(ImageInterface<Float>& image,
160 : ImageBeamSet& beams, Bool normalizeVolume){
161 :
162 0 : Int nbeams=beams.shape()(0);
163 :
164 0 : Int freqAx=CoordinateUtil::findSpectralAxis(image.coordinates());
165 0 : Int nchan=image.shape()(freqAx);
166 0 : if((nchan != nbeams) || (nchan==1)){
167 0 : GaussianBeam elbeam=beams(0, 0);
168 0 : Convolve(image, elbeam, normalizeVolume);
169 0 : }
170 : else{
171 0 : IPosition blc=image.shape();
172 0 : blc=0;
173 0 : IPosition trc=image.shape()-1;
174 0 : for (Int k=0; k < nchan; ++k){
175 0 : blc[freqAx]=k;
176 0 : trc[freqAx]=k;
177 0 : Slicer slc(blc, trc, Slicer::endIsLast);
178 0 : SubImage<Float> subim(image, slc, true);
179 0 : GaussianBeam elbeam=beams(k,0);
180 0 : Convolve(subim, elbeam, normalizeVolume);
181 :
182 0 : }
183 0 : }
184 0 : }
185 :
186 3130 : void StokesImageUtil::Convolve(ImageInterface<Float>& image,
187 : GaussianBeam& beam,
188 : Bool normalizeVolume)
189 : {
190 3130 : Convolve(image, beam.getMajor().get("arcsec").getValue(),
191 6260 : beam.getMinor().get("arcsec").getValue(),
192 6260 : beam.getPA().get("deg").getValue(),
193 : normalizeVolume);
194 3130 : }
195 :
196 3130 : void StokesImageUtil::Convolve(ImageInterface<Float>& image, Float bmaj,
197 : Float bmin, Float bpa,
198 : Bool normalizeVolume) {
199 :
200 3130 : Vector<Float> beam(3);
201 3130 : beam(0)=bmaj*C::arcsec;
202 3130 : beam(1)=bmin*C::arcsec;
203 3130 : beam(2)=(bpa+90.0)*C::degree;
204 : PagedImage<Float>
205 6260 : cleanpsf(IPosition(4, image.shape()(0), image.shape()(1), 1, 1),
206 : image.coordinates(),
207 12520 : File::newUniqueName("./","imagesolver::cleanpsf").originalName());
208 3130 : cleanpsf.table().markForDelete();
209 3130 : MakeGaussianPSF(cleanpsf, beam, normalizeVolume);
210 3130 : Convolve(image, cleanpsf);
211 3130 : }
212 :
213 0 : void StokesImageUtil::MakeGaussianPSF(ImageInterface<Float>& psf,
214 : Quantity& bmaj, Quantity& bmin,
215 : Quantity& bpa,
216 : Bool normalizeVolume)
217 : {
218 0 : Vector<Float> beam(3);
219 0 : beam(0)=bmaj.get("arcsec").getValue();
220 0 : beam(1)=bmin.get("arcsec").getValue();
221 0 : beam(2)=bpa.get("deg").getValue();
222 0 : MakeGaussianPSF(psf, beam, normalizeVolume);
223 0 : }
224 :
225 : // Make an image into a Gaussian PSF
226 5774 : void StokesImageUtil::MakeGaussianPSF(ImageInterface<Float>& psf, Vector<Float>& beam,
227 : Bool normalizeVolume) {
228 :
229 5774 : Int nx=psf.shape()(0);
230 5774 : Int ny=psf.shape()(1);
231 5774 : Matrix<Float> ipsf(nx, ny);
232 :
233 5774 : Double cospa=cos(beam(2));
234 5774 : Double sinpa=sin(beam(2));
235 5774 : Vector<Double> rp=psf.coordinates().referencePixel();
236 5774 : Vector<Double> d=psf.coordinates().increment();
237 5774 : Double refi=rp(0);
238 5774 : Double refj=rp(1);
239 5774 : AlwaysAssert(beam(0)>0.0,AipsError);
240 5774 : AlwaysAssert(beam(1)>0.0,AipsError);
241 : Double sbmaj, sbmin;
242 :
243 : // ---old
244 : // Assumes that the cell sizes are the same in both directions!
245 : //sbmaj=4.0*log(2.0)*square(d(0)/beam(0));
246 : //sbmin=4.0*log(2.0)*square(d(0)/beam(1));
247 : // ---old
248 :
249 : // Evaluate the Gaussian beam using real-world coordinates, to handle rectangular pixels (cas-7171)
250 5774 : sbmaj=4.0*log(2.0)*square(1.0/beam(0));
251 5774 : sbmin=4.0*log(2.0)*square(1.0/beam(1));
252 :
253 5774 : Vector<Double>fd(fabs(d));
254 :
255 5774 : Float volume=0.0;
256 1155612 : for (Int j=0;j<ny;j++) {
257 655406474 : for (Int i=0;i<nx;i++) {
258 : // ---old
259 : // Double x = cospa * (Double(i)-refi) + sinpa * (Double(j)-refj);
260 : // Double y = - sinpa * (Double(i)-refi) + cospa * (Double(j)-refj);
261 : // ---old
262 654256636 : Double x = cospa * (Double(i)-refi)*fd(0) + sinpa * (Double(j)-refj)*fd(1);
263 654256636 : Double y = - sinpa * (Double(i)-refi)*fd(0) + cospa * (Double(j)-refj)*fd(1);
264 654256636 : Double radius = sbmaj*square(x) + sbmin*square(y);
265 654256636 : if (radius<20.) {
266 24414267 : ipsf(i,j) = exp(-radius);
267 24414267 : volume+=ipsf(i,j);
268 : }
269 : else {
270 629842369 : ipsf(i,j)=0.;
271 : }
272 : }
273 : }
274 5774 : if(normalizeVolume) ipsf/=volume;
275 5774 : psf.putSlice(ipsf, IPosition(psf.ndim(),0), IPosition(psf.ndim(),1));
276 :
277 5774 : }
278 :
279 2 : Float StokesImageUtil::normalizePSF(casacore::ImageInterface<casacore::Float>& psf){
280 2 : AlwaysAssert(psf.ndim()==4,AipsError);
281 :
282 2 : Float retval=-1e10;
283 2 : Vector<Int> map;
284 2 : AlwaysAssert(StokesMap(map, psf.coordinates()), AipsError);
285 2 : Int nx = psf.shape()(map(0));
286 2 : Int ny = psf.shape()(map(1));
287 2 : if(nx < 32 || ny < 32)
288 0 : throw(AipsError("Will not normalize images less than 32 pixels"));
289 2 : Int npol = psf.shape()(map(2));
290 2 : Int nchan = psf.shape()(map(3));
291 2 : IPosition blc(4);
292 2 : IPosition trc(4);
293 2 : IPosition blc8(4,0,0,0,0);
294 2 : IPosition trc8(4,0,0,0,0);
295 2 : blc(map(0))=0; blc(map(1))=0;
296 2 : trc(map(0))=nx-1; trc(map(1))=ny-1;
297 2 : blc8(map(0))=nx/2-nx/16; blc8(map(1))=ny/2-ny/16;
298 2 : trc8(map(0))=nx/2+nx/16; trc8(map(1))=ny/2+ny/16;
299 2 : Slicer inner8(blc8, trc8, Slicer::endIsLast);
300 6 : for (Int k=0; k < nchan ; ++k){
301 4 : blc(map(3))=k;
302 4 : trc(map(3))=k;
303 8 : for (Int j=0; j < npol; ++j){
304 4 : blc(map(2))=j;
305 4 : trc(map(2))=j;
306 4 : Slicer sl(blc, trc, Slicer::endIsLast);
307 4 : SubImage<Float> psfSub(psf, sl, True);
308 4 : Float maxCenter=max(psfSub.getSlice(inner8));
309 4 : if(maxCenter > retval)
310 2 : retval=maxCenter;
311 4 : if(maxCenter !=0){
312 4 : psfSub.copyData( (LatticeExpr<Float>)(psfSub/maxCenter));
313 : }
314 : else{
315 0 : psfSub.set(0.0);
316 : }
317 4 : }
318 : }
319 2 : return retval;
320 2 : }
321 :
322 :
323 : // Zero specified elements of a Stokes image
324 0 : void StokesImageUtil::Zero(ImageInterface<Float>& image, Vector<Bool>& mask) {
325 :
326 0 : AlwaysAssert(image.ndim()==4,AipsError);
327 0 : AlwaysAssert(mask.shape()(0)==4,AipsError);
328 :
329 0 : Vector<Int> map;
330 0 : AlwaysAssert(StokesMap(map, image.coordinates()), AipsError);
331 :
332 0 : Int nx = image.shape()(map(0));
333 0 : Int ny = image.shape()(map(1));
334 0 : Int npol = image.shape()(map(2));
335 0 : Int nchan = image.shape()(map(3));
336 :
337 0 : LatticeStepper ls(image.shape(), IPosition(4, nx, 1, 1, 1),
338 0 : IPosition(4, map(0), map(1), map(2), map(3)));
339 0 : LatticeIterator<Float> li(image, ls);
340 :
341 : // Loop over all planes
342 0 : li.reset();
343 0 : for (Int chan=0;chan<nchan;chan++) {
344 0 : for(Int pol=0;pol<npol;pol++) {
345 0 : for(Int iy=0;iy<ny;iy++,li++) {
346 0 : if(mask(pol)) {
347 0 : li.woVectorCursor()=0.0;
348 : }
349 : }
350 : }
351 : }
352 :
353 0 : }
354 :
355 0 : void StokesImageUtil::MaskFrom(ImageInterface<Float>& mask,
356 : ImageInterface<Float>& image,
357 : const Quantity& threshold)
358 : {
359 0 : MaskFrom(mask, image, threshold.get("Jy").getValue());
360 0 : }
361 :
362 0 : void StokesImageUtil::MaskFrom(ImageInterface<Float>& mask,
363 : ImageInterface<Float>& image,
364 : const Double thres)
365 : {
366 0 : AlwaysAssert(image.ndim()==4,AipsError);
367 0 : AlwaysAssert(mask.shape()==image.shape(),AipsError);
368 :
369 0 : Vector<Int> map;
370 0 : AlwaysAssert(StokesMap(map, image.coordinates()), AipsError);
371 :
372 0 : Int nx = image.shape()(map(0));
373 0 : Int npol = image.shape()(map(2));
374 :
375 0 : LatticeStepper ls(image.shape(), IPosition(4, nx, 1, npol, 1),
376 0 : IPosition(4, map(0), map(1), map(2), map(3)));
377 0 : RO_LatticeIterator<Float> li(image, ls);
378 0 : LatticeIterator<Float> mli(mask, ls);
379 :
380 : // Loop over all planes
381 0 : for (li.reset(),mli.reset();!li.atEnd()&&!mli.atEnd();li++,mli++) {
382 0 : if(npol>1) {
383 0 : mli.rwMatrixCursor()=0.0;
384 0 : for (Int ix=0;ix<nx;ix++) {
385 0 : if(li.matrixCursor()(ix,0)>thres) {
386 0 : for (Int pol=0;pol<npol;pol++) mli.rwMatrixCursor()(ix,pol)=1.0;
387 : }
388 : }
389 : }
390 : else {
391 0 : mli.rwVectorCursor()=0.0;
392 0 : for (Int ix=0;ix<nx;ix++) {
393 0 : if(li.vectorCursor()(ix)>thres) {
394 0 : mli.rwVectorCursor()(ix)=1.0;
395 : }
396 : }
397 : }
398 : }
399 0 : }
400 :
401 0 : void StokesImageUtil::MaskOnStokesI(ImageInterface<Float>& image,
402 : const Quantity& threshold)
403 : {
404 :
405 0 : Double thres=threshold.get("Jy").getValue();
406 :
407 0 : AlwaysAssert(image.ndim()==4,AipsError);
408 :
409 0 : Vector<Int> map;
410 0 : AlwaysAssert(StokesMap(map, image.coordinates()), AipsError);
411 :
412 0 : Int nx = image.shape()(map(0));
413 0 : Int npol = image.shape()(map(2));
414 :
415 0 : LatticeStepper ls(image.shape(), IPosition(4, nx, 1, npol, 1),
416 0 : IPosition(4, map(0), map(1), map(2), map(3)));
417 0 : LatticeIterator<Float> li(image, ls);
418 :
419 : // Loop over all planes
420 0 : for (li.reset();!li.atEnd();li++) {
421 0 : for (Int ix=0;ix<nx;ix++) {
422 0 : if(npol>1) {
423 0 : if(li.matrixCursor()(ix,0)<thres) {
424 0 : for (Int pol=0;pol<npol;pol++) li.rwMatrixCursor()(ix,pol)=0.0;
425 : }
426 : }
427 : else {
428 0 : if(li.vectorCursor()(ix)<thres) {
429 0 : li.rwVectorCursor()(ix)=0.0;
430 : }
431 : }
432 : }
433 : }
434 :
435 0 : }
436 :
437 : // Constrain (a) Stokes I to be > 0 and (b) I > abs(V) and
438 : // (c) fractional pol<1.0
439 0 : void StokesImageUtil::Constrain(ImageInterface<Float>& image) {
440 :
441 0 : AlwaysAssert(image.ndim()==4,AipsError);
442 :
443 0 : Vector<Int> map;
444 0 : AlwaysAssert(StokesMap(map, image.coordinates()), AipsError);
445 :
446 0 : Int nx = image.shape()(map(0));
447 0 : Int npol = image.shape()(map(2));
448 :
449 0 : LatticeStepper ls(image.shape(), IPosition(4, nx, 1, npol, 1),
450 0 : IPosition(4, map(0), map(2), map(1), map(3)));
451 0 : LatticeIterator<Float> li(image, ls);
452 :
453 : // Loop over all planes
454 0 : for (li.reset();!li.atEnd();li++) {
455 :
456 : // Make a copy of the cursor
457 0 : Matrix<Float>& limage=li.rwMatrixCursor();
458 :
459 0 : for (Int ix=0;ix<nx;ix++) {
460 : // I >=0
461 0 : Bool zero=false;
462 0 : if(limage(ix,0)<0.0) {
463 0 : zero=true;
464 : }
465 : // I >= abs(V)
466 0 : else if (npol==2) {
467 0 : if(limage(ix,0)<abs(limage(ix,1))) zero=true;
468 : }
469 : // fractional pol < 1.0
470 0 : else if (npol==4) {
471 0 : StokesVector sv(limage(ix,0),limage(ix,1),
472 0 : limage(ix,2),limage(ix,3));
473 0 : if(sv.minEigenValue()<0.0) zero=true;
474 : }
475 0 : if(zero) for (Int pol=0;pol<npol;pol++) limage(ix,pol)=0.0;
476 : }
477 : }
478 0 : }
479 :
480 : //See CAS-13022 for FitGaussianPSF algorithm details
481 42 : Bool StokesImageUtil::FitGaussianPSF(ImageInterface<Float>& psf, ImageBeamSet& elbeam, Float psfcutoff){
482 :
483 42 : Bool retval=true;
484 42 : Int freqAx=CoordinateUtil::findSpectralAxis(psf.coordinates());
485 42 : Vector<Stokes::StokesTypes> whichPols;
486 42 : Int polAx=CoordinateUtil::findStokesAxis(whichPols, psf.coordinates());
487 42 : Int nchan=psf.shape()(freqAx);
488 42 : IPosition blc=psf.shape();
489 42 : blc=0;
490 42 : IPosition trc=psf.shape()-1;
491 42 : trc[polAx]=0;
492 42 : elbeam=ImageBeamSet(nchan, 1);
493 42 : IPosition ipos(2,0,0);
494 42 : Matrix<GaussianBeam> tempBeam(nchan,1);
495 42 : Vector<Bool> fitted(nchan, false);
496 :
497 84 : for (Int k=0; k < nchan;++k){
498 42 : blc[freqAx]=k;
499 42 : trc[freqAx]=k;
500 42 : Slicer slc(blc, trc, Slicer::endIsLast);
501 42 : SubImage<Float> subpsf(psf, slc, false);
502 : try{
503 42 : fitted(k)=FitGaussianPSF(subpsf, tempBeam(k,0), psfcutoff);
504 : }
505 0 : catch (AipsError x_error){
506 0 : Int ik=k;
507 0 : fitted(k)=false;
508 0 : while ((ik > 0) && !fitted(k)){
509 0 : if(fitted(ik-1)){
510 0 : fitted(k)=true;
511 0 : tempBeam(k,0)=tempBeam(ik-1, 0);
512 : }
513 0 : --ik;
514 : }
515 0 : }
516 42 : ipos(0)=k;
517 42 : }
518 42 : Float maxMaj=0.0;
519 42 : Float minMaj=C::flt_max;
520 42 : Int posMax=0;
521 42 : DirectionCoordinate directionCoord=psf.coordinates().directionCoordinate(psf.coordinates().findCoordinate(Coordinate::DIRECTION));
522 42 : String dirunit=directionCoord.worldAxisUnits()(0);
523 84 : for(Int k=0; k < nchan; ++k){
524 42 : if(fitted(k) && (maxMaj < tempBeam(k,0).getMajor(dirunit))){
525 42 : maxMaj=tempBeam(k,0).getMajor(dirunit);
526 42 : posMax=k;
527 : }
528 42 : if(fitted(k) && (minMaj > tempBeam(k,0).getMajor(dirunit))){
529 42 : minMaj=tempBeam(k,0).getMajor(dirunit);
530 42 : posMax=k;
531 : }
532 : }
533 : //If the beams are within 0.5 pixel the same resolution then
534 : //who cares to have a beam per plane
535 42 : if(abs(minMaj-maxMaj) < 0.5*abs(directionCoord.increment()(0))){
536 42 : GaussianBeam theBeam=tempBeam(posMax,0);
537 42 : tempBeam.resize(1,1);
538 42 : tempBeam(0,0)=theBeam;
539 42 : fitted.resize(1);
540 42 : fitted=true;
541 42 : }
542 42 : if(!anyTrue(fitted)){
543 0 : retval=false;
544 0 : throw(AipsError("No valid fits were found to PSF"));
545 : }
546 42 : if(!allTrue(fitted)){
547 0 : for (Int k=0; k < nchan;++k){
548 0 : int ik=k;
549 0 : while((ik < (nchan-1)) && !fitted(k)){
550 0 : if(fitted(ik+1)){
551 0 : fitted(k)=true;
552 0 : tempBeam(k,0)=tempBeam(ik+1, 0);
553 : }
554 0 : ++ik;
555 : }
556 : }
557 : }
558 42 : elbeam=ImageBeamSet(tempBeam.shape());
559 42 : elbeam.setBeams(tempBeam);
560 :
561 42 : return retval;
562 42 : }
563 :
564 2940 : Bool StokesImageUtil::FitGaussianPSF(ImageInterface<Float>& psf,
565 : GaussianBeam& beam, Float psfcutoff)
566 : {
567 2940 : Vector<Float> vbeam(3, 0.0);
568 2940 : Bool status=true;
569 2940 : if(!FitGaussianPSF(psf, vbeam, psfcutoff)) status=false;
570 5880 : beam = GaussianBeam(
571 5880 : Quantity(abs(vbeam[0]),"arcsec"),
572 5880 : Quantity(abs(vbeam[1]),"arcsec"),
573 5880 : Quantity(vbeam[2],"deg")
574 2940 : );
575 2940 : return status;
576 :
577 2940 : }
578 :
579 :
580 5880 : void StokesImageUtil::FindNpoints(Int& npoints, IPosition& blc, IPosition& trc, Int nrow, Float amin, Int px, Int py, Vector<Double>& deltas, Matrix<Double>& x , Vector<Double>& y, Vector<Double>& sigma, Matrix<Float>& lpsf){
581 :
582 5880 : Int maxnpoints=(2*nrow+1)*(2*nrow+1);
583 5880 : Matrix<Double> ix(maxnpoints,2);
584 5880 : Vector<Double> iy(maxnpoints), isigma(maxnpoints);
585 5880 : ix=0.0; iy=0.0; isigma=0.0;
586 5880 : npoints = 0;
587 :
588 : //IPosition blc(2), trc(2);
589 : //blc = 0; trc = 0;
590 :
591 5880 : blc(0) = lpsf.shape()(0)-1;
592 5880 : blc(1) = lpsf.shape()(1)-1;
593 5880 : trc(0) = 0;
594 5880 : trc(1) = 0;
595 :
596 5880 : IPosition psf_shape = lpsf.shape();
597 :
598 : //we sample the central part of a, 2*nrow+1 by 2*nrow+1
599 :
600 5880 : Int iflip = 1;
601 5880 : Int jflip = 1;
602 : // loop through rows. Include both above and below in case
603 : // we are fitting an image feature
604 :
605 17595 : for(Int jlo = 0;jlo<2;jlo++) {
606 11760 : jflip*=-1;
607 : // loop from 0 to nrow then from 1 to nrow
608 261077 : for(Int j = jlo;j<=nrow;j++) {
609 : // the current row under consideration
610 249362 : Int jrow = py + j*jflip;
611 : // loop down row doing alternate halves. work our way
612 : // out from the center until we cross threshold
613 : // don't include any sidelobes!
614 740855 : for(Int ilo=0;ilo<2;ilo++) {
615 495131 : iflip*=-1;
616 : // start at center row this may or may not be in the lobe,
617 : // if it's narrow and the pa is near 45 degrees
618 :
619 :
620 495131 : if((jrow > (psf_shape(1)-1)) || (jrow < 0 ) ) break;
621 : //cout << "11.** " << maxnpoints << ",*,"<< npoints << ",*," << px << ",*," << jrow << endl;
622 491538 : Bool inlobe = lpsf(px,jrow)>amin;
623 :
624 : //cout << "12.** " << maxnpoints << ",*,"<< npoints << ",*," << px << ",*," << jrow << endl;
625 8493421 : for(Int i = ilo;i<=nrow;i++) {
626 8242763 : if(npoints < maxnpoints){
627 8242763 : Int irow = px + i*iflip;
628 : // did we step out of the lobe?
629 :
630 8242763 : if((irow > (psf_shape(0)-1)) || (irow < 0 ) ) break;
631 : //cout << "irow , jrow " << irow << ",*," << jrow << ",*," << lpsf.shape() << ",*," << px << ",*," << py << endl;
632 8173845 : if (inlobe&&(lpsf(irow,jrow)<amin)) break;
633 8001928 : if (lpsf(irow,jrow)>amin) {
634 : //cout << "$%$irow , jrow " << lpsf.shape() << ",*,"<< irow << "," << jrow << ",*, "<< irow - px << ",*," << jrow - py << ",*, " << lpsf(irow,jrow) << ",*, " << inlobe << ",*, " << npoints << endl;
635 1287926 : inlobe = true;
636 : // the sign on the ra can cause problems. we just fit
637 : // for what the beam "looks" like here, and worry about
638 : // it later.
639 1287926 : ix(npoints,0) = (irow-px)*abs(deltas(0));
640 1287926 : ix(npoints,1) = (jrow-py)*abs(deltas(1));
641 1287926 : iy(npoints) = lpsf(irow,jrow);
642 : //cout << "1.**" << endl;
643 :
644 1287926 : if(blc(0) > irow) blc(0) = irow;
645 1287926 : if(blc(1) > jrow) blc(1) = jrow;
646 :
647 1287926 : if(trc(0) < irow) trc(0) = irow;
648 1287926 : if(trc(1) < jrow) trc(1) = jrow;
649 : //cout << "2.**" << endl;
650 :
651 1287926 : isigma(npoints) = 1.0;
652 1287926 : ++npoints;
653 : // cout << "3.** " << maxnpoints << ",*,"<< npoints << endl;
654 1287926 : if(npoints > (maxnpoints-1)) {
655 45 : inlobe=false;
656 : //cout << "Over max points .** " << maxnpoints << ",*,"<< npoints << endl;
657 : //break;
658 45 : goto endSearch;
659 : }
660 : // cout << "5.** " << maxnpoints << ",*,"<< npoints << endl;
661 : }
662 : }
663 : }
664 : }
665 : }
666 : }
667 :
668 5835 : endSearch:
669 :
670 : //cout << "2. blc, trc " << blc << ",*," << trc << endl;
671 : //Vector<Double> y(npoints), sigma(npoints);
672 : //Matrix<Double> x(npoints,2);
673 : //cout << "hallo " << endl;
674 :
675 5880 : y.resize(npoints);
676 5880 : x.resize(npoints,2);
677 5880 : sigma.resize(npoints);
678 :
679 : //cout << "2. findNPoints " << npoints << " " << amin << endl;
680 :
681 1293806 : for (Int ip=0;ip<npoints;ip++) {
682 1287926 : x(ip,0)=ix(ip,0); x(ip,1)=ix(ip,1);
683 1287926 : y(ip)=iy(ip);
684 1287926 : sigma(ip)=isigma(ip);
685 1287926 : if(!(isigma(ip)>0.0)) break;
686 : }
687 :
688 : //Ensure that it is square
689 5880 : if(blc(0) > blc(1)){
690 1112 : blc(0) = blc(1);
691 : }else{
692 4768 : blc(1) = blc(0);
693 : }
694 :
695 5880 : if(trc(0) > trc(1)){
696 2431 : trc(1) = trc(0);
697 : }else{
698 3449 : trc(0) = trc(1);
699 : }
700 :
701 :
702 5880 : }
703 :
704 2940 : void StokesImageUtil::ResamplePSF(Matrix<Float>& psf, Int& oversampling, Matrix<Float>& resampledPsf, String& InterpMethod)
705 : {
706 2940 : Int nx = psf.shape()(0);
707 2940 : Int ny = psf.shape()(1);
708 2940 : Int nxRe = nx*oversampling - oversampling + 1;
709 2940 : Int nyRe = ny*oversampling - oversampling + 1;
710 :
711 2940 : Vector<Double> pos(2);
712 2940 : resampledPsf.resize(nxRe,nyRe);
713 :
714 : Interpolate2D::Method method;
715 2940 : if(InterpMethod == "NEAREST") method = Interpolate2D::NEAREST;
716 2940 : if(InterpMethod == "LINEAR") method = Interpolate2D::LINEAR;
717 2940 : if(InterpMethod == "CUBIC") method = Interpolate2D::CUBIC;
718 2940 : if(InterpMethod == "LANCZOS") method = Interpolate2D::LANCZOS;
719 :
720 2940 : Interpolate2D resampleInterp(method);
721 : //Interpolate2D resampleInterp(Interpolate2D::CUBIC);
722 :
723 : Bool ok;
724 : Float result;
725 :
726 129526 : for (Int i=0; i < nxRe ; ++i){
727 5977708 : for (Int j=0; j < nyRe ; ++j){
728 5851122 : pos(0) = (Float) i/(Float) oversampling;
729 5851122 : pos(1) = (Float) j/(Float) oversampling;
730 :
731 5851122 : ok = resampleInterp.interp(result, pos, psf);
732 : //cout << "pos is" << pos << " result " << result << " psf " << psf(i,j)<< endl;
733 5851122 : resampledPsf(i,j) = result;
734 : }
735 : }
736 :
737 :
738 2940 : }
739 :
740 :
741 2940 : Bool StokesImageUtil::FitGaussianPSF(ImageInterface<Float>& psf, Vector<Float>& beam, Float psfcutoff) {
742 :
743 5880 : LogIO os(LogOrigin("StokesImageUtil", "FitGaussianPSF()",WHERE));
744 2940 : os << LogIO::DEBUG1 << "Psfcutoff is " << psfcutoff << LogIO::POST;
745 :
746 2940 : Float npix = 20;
747 2940 : Int expand_pixel = 5;
748 2940 : Int target_npoints = 3001;
749 2940 : String InterpMethod = "CUBIC";
750 :
751 : //##########################################################
752 :
753 2940 : Vector<Double> deltas;
754 2940 : deltas=psf.coordinates().increment();
755 :
756 2940 : Vector<Int> map;
757 2940 : AlwaysAssert(StokesMap(map, psf.coordinates()), AipsError);
758 :
759 2940 : Int nx = psf.shape()(map(0));
760 2940 : Int ny = psf.shape()(map(1));
761 : //// Int nchan = psf.shape()(map(3)); // apparently not used and no side effects
762 2940 : Int px=0;
763 2940 : Int py=0;
764 2940 : Float bamp=0;
765 2940 : Matrix<Float> lpsf;
766 2940 : locatePeakPSF(psf, px, py, bamp, lpsf);
767 :
768 : //check if peak is ZERO
769 2940 : if( bamp < 1e-07 ) {
770 0 : throw(AipsError("Psf peak is zero"));
771 : }
772 :
773 : //check if peak is outside inner quarter
774 2940 : if(px < nx/4.0 || px > 3.0*nx/4.0 || py < ny/4.0 || py > 3.0*ny/4.0) {
775 0 : throw(AipsError("Peak of psf is outside the inner quarter of defined image"));
776 : }
777 :
778 :
779 2940 : if((bamp > 1.1) || (bamp < 0.9)) // No warning if 10% error in PSF peak
780 0 : os << LogIO::WARN << "Peak of PSF is " << bamp << LogIO::POST;
781 :
782 :
783 2940 : static const Float fdiam = 2.5*abs(deltas(0))/C::arcsec;
784 : try{
785 :
786 :
787 2940 : if(bamp==0.0) {
788 0 : beam[0] = fdiam;
789 0 : beam[1] = fdiam;
790 0 : beam[2] = 0.0;
791 0 : os << LogIO::WARN << "Could not find peak " << LogIO::POST;
792 0 : return false;
793 : }
794 :
795 :
796 :
797 2940 : lpsf/=bamp; //Normalize
798 :
799 : // The selection code attempts to avoid including any sidelobes, even
800 : // if they exceed the threshold, by starting from the center column and
801 : // working out, exiting the loop when it crosses the threshold. It
802 : // assumes that the first time it finds a "good" point starting from
803 : // the center and working out that it's in the main lobe. Narrow,
804 : // sharply ringing beams inclined at 45 degrees will confuse it, but
805 : // that's even more pathological than most VLBI beams.
806 :
807 2940 : Float amin=psfcutoff;
808 2940 : Int nrow=npix;
809 :
810 : //we sample the central part of a, 2*nrow+1 by 2*nrow+1
811 2940 : Bool converg=false;
812 2940 : Vector<Double> solution;
813 2940 : Int kounter=0;
814 :
815 5880 : while(amin > 0.009 && !converg && (kounter < 50)){
816 2940 : kounter+=1;
817 2940 : Int npoints=0;
818 2940 : Vector<Double> y, sigma;
819 2940 : Matrix<Double> x;
820 2940 : IPosition blc(2), trc(2);
821 :
822 2940 : FindNpoints(npoints, blc, trc, nrow, amin, px, py, deltas, x , y, sigma, lpsf);
823 2940 : os << LogIO::DEBUG1 << "First FindNpoints is " << npoints << LogIO::POST;
824 :
825 2940 : blc = blc-expand_pixel;
826 2940 : trc = trc+expand_pixel;
827 :
828 2940 : if(blc(0) < 0) blc(0)=0;
829 2940 : if(blc(1) < 0) blc(1)=0;
830 2940 : if(trc(0) >= lpsf.shape()(0)) trc(0)=lpsf.shape()(0)-1;
831 2940 : if(trc(1) >= lpsf.shape()(1)) trc(1)=lpsf.shape()(1)-1;
832 :
833 2940 : Matrix<Float> lpsfWindowed = lpsf(blc,trc);
834 2940 : os << LogIO::DEBUG1 << "The windowed Psf shape is " << lpsfWindowed.shape() << LogIO::POST;
835 2940 : Matrix<Float> resampledPsf;
836 :
837 2940 : Int oversampling = (Int) sqrt(target_npoints/(lpsfWindowed.shape()(0)*lpsfWindowed.shape()(1)));
838 2940 : if(oversampling == 0){
839 9 : oversampling = 1;
840 : }
841 2940 : os << LogIO::DEBUG1 << "The oversampling is " << oversampling << LogIO::POST;
842 :
843 2940 : ResamplePSF(lpsfWindowed, oversampling, resampledPsf,InterpMethod);
844 2940 : os << LogIO::DEBUG1 << "The resampled windowed Psf shape is " << lpsfWindowed.shape() << LogIO::POST;
845 :
846 : Float minVal, maxVal;
847 2940 : IPosition minPos(2);
848 2940 : IPosition maxPos(2);
849 2940 : minMax(minVal, maxVal, minPos, maxPos, resampledPsf);
850 2940 : resampledPsf = resampledPsf/maxVal;
851 :
852 2940 : Vector<Double> resampledDeltas = deltas/ (Double) oversampling;
853 :
854 : Int minLen;
855 2940 : if((trc(0) - blc(0)) > (trc(1) - blc(1)) ){
856 0 : minLen = trc(1) - blc(1) + 1;
857 : }else{
858 2940 : minLen = trc(0) - blc(0) + 1;
859 : }
860 :
861 2940 : Int nrowRe = (Int) (oversampling*minLen - 1)/2;
862 2940 : FindNpoints(npoints, blc, trc, nrowRe, amin, maxPos(0), maxPos(1), resampledDeltas, x , y, sigma, resampledPsf);
863 2940 : os << LogIO::DEBUG1 << "Second FindNpoints is " << npoints << LogIO::POST;
864 :
865 2940 : Gaussian2D<AutoDiff<Double> > gauss2d;
866 2940 : gauss2d[0] = 1.0; //Height of Gaussian
867 2940 : gauss2d[1] = 0.0; //The center of the Gaussian in the x direction
868 2940 : gauss2d[2] = 0.0; //The center of the Gaussian in the y direction.
869 : //gauss2d[3] = 2.5*abs(resampledDeltas(0)); //The width (FWHM) of the Gaussian on one axis.
870 2940 : gauss2d[3] = 2.5*abs(deltas(0));
871 2940 : gauss2d[4] = 0.5; //A modified axial ratio.
872 :
873 : // Fix height and center
874 2940 : gauss2d.mask(0) = false;
875 2940 : gauss2d.mask(1) = false;
876 2940 : gauss2d.mask(2) = false;
877 :
878 : // CAS-13515: Fitting sometimes fails with "NonLinearFitLM: error in loop solution". This occurs in casacore/scimath/Fitting/NonLinearFitLM.tcc LSQFit::invertRect() due to a matrix that can not be inverted.
879 : // This problem is solved by retrying fit with a different position angle.
880 2940 : Bool loopSolutionFound=false;
881 2940 : Int retryCounter = 0;
882 2940 : Double posAng = 1.0; //Initial position angle.
883 2940 : Int nRetries = 10;
884 5880 : while(!loopSolutionFound && retryCounter < nRetries){
885 : //Create casacore fitter and fit Gaussian to subset of PSF interpolated values.
886 2940 : gauss2d[5] = posAng; //The position angle.
887 2940 : NonLinearFitLM<Double> fitter;
888 : // Set maximum number of iterations to 1000
889 2940 : fitter.setMaxIter(1000);
890 :
891 : // Set converge criteria. Default is 0.001
892 2940 : fitter.setCriteria(0.0001);
893 :
894 : // Set the function and initial values
895 2940 : fitter.setFunction(gauss2d);
896 :
897 : try{
898 : // The current parameter values are used as the initial guess.
899 2940 : solution = fitter.fit(x, y, sigma);
900 2940 : loopSolutionFound = true;
901 0 : }catch(AipsError x_error){
902 0 : loopSolutionFound = false;
903 0 : retryCounter++;
904 0 : posAng = 1.0 + retryCounter*M_PI/nRetries;
905 0 : os << LogIO::NORMAL3 << "Fit failed, another atempt will be made with position angle " << posAng << " rad." << LogIO::POST;
906 0 : }
907 2940 : converg=fitter.converged();
908 2940 : }
909 :
910 2940 : if(!loopSolutionFound)
911 : {
912 0 : throw(AipsError(String("Error in StokesImageUtil::FitGaussianPSF: error in loop solution.")));
913 : }
914 :
915 2940 : if (!converg) {
916 0 : beam(0)=2.5*abs(deltas(0))/C::arcsec;
917 0 : beam(1)=2.5*abs(deltas(0))/C::arcsec;
918 0 : beam(2)=0.0;
919 :
920 : //fit did not coverge...reduce the minimum i.e expand the field a bit
921 0 : amin/=1.5;
922 :
923 0 : os << LogIO::WARN << "The fit did not coverge; another atempt will be made with psfcutoff " << amin << LogIO::POST;
924 : }
925 :
926 2940 : }
927 2940 : if (converg) {
928 2940 : if (abs(solution(4))>1.0) {
929 141 : beam(0)=abs(solution(3)*solution(4))/C::arcsec;
930 141 : beam(1)=abs(solution(3))/C::arcsec;
931 141 : beam(2)=solution(5)/C::degree-90.0;
932 : } else {
933 2799 : beam(0)=abs(solution(3))/C::arcsec;
934 2799 : beam(1)=abs(solution(3)*solution(4))/C::arcsec;
935 2799 : beam(2)=solution(5)/C::degree;
936 : }
937 :
938 2940 : beam(2)=fmod(beam(2), Float(360.0));
939 : //while (abs(beam(2)/180.0)> 1) {
940 : // if (beam(2) > 180.0) beam(2)-=360.0;
941 : // else beam(2)+=360.0;
942 : //}
943 : //CAS-8627: adjust pa to the range -90d, +90d
944 4799 : while (abs(beam(2)/90.0)> 1) {
945 1859 : if (beam(2) > 270.0) beam(2)-=360.0;
946 1840 : else if (beam(2) > 90.0) beam(2) -=180.0;
947 35 : else if (beam(2) < -270.0) beam(2) +=360.0;
948 31 : else beam(2)+=180.0;
949 : }
950 :
951 2940 : return true;
952 : }
953 : else os << LogIO::WARN << "The fit did not coverge; check your PSF" <<
954 0 : LogIO::POST;
955 0 : return false;
956 :
957 2940 : } catch (AipsError x_error) {
958 0 : beam[0] = fdiam;
959 0 : beam[1] = fdiam;
960 0 : beam[2] = 0.0;
961 : os << LogIO::SEVERE << "Caught Exception: "<< x_error.getMesg() <<
962 0 : LogIO::POST;
963 0 : return false;
964 0 : }
965 :
966 2940 : }
967 :
968 2940 : void StokesImageUtil::locatePeakPSF(ImageInterface<Float>& in, Int& xpos, Int& ypos, Float& amp, Matrix<Float>& lpsf){
969 2940 : Vector<Int> map;
970 2940 : AlwaysAssert(StokesMap(map, in.coordinates()), AipsError);
971 :
972 2940 : Int nx = in.shape()(map(0));
973 2940 : Int ny = in.shape()(map(1));
974 2940 : Int nchan = in.shape()(map(3));
975 2940 : xpos=0;
976 2940 : ypos=0;
977 : Float psfmax;
978 :
979 2940 : amp=0;
980 :
981 :
982 2940 : IPosition oneplane(in.ndim(), nx, ny, 1, 1);
983 5880 : LatticeStepper psfls(in.shape(), oneplane,
984 5880 : IPosition(4,map(0),map(1),map(3),map(2)));
985 2940 : RO_LatticeIterator<Float> psfli(in,psfls);
986 5880 : while (nchan >= 1 && amp < 0.9){
987 2940 : lpsf=psfli.matrixCursor();
988 :
989 2940 : IPosition psfposmax(lpsf.ndim());
990 2940 : IPosition psfposmin(lpsf.ndim());
991 : Float psfmin;
992 :
993 2940 : minMax(psfmin, psfmax, psfposmin, psfposmax, lpsf);
994 :
995 2940 : xpos=psfposmax(0);
996 2940 : ypos=psfposmax(1);
997 :
998 :
999 :
1000 2940 : amp=lpsf(xpos,ypos);
1001 :
1002 2940 : ++psfli;
1003 2940 : --nchan;
1004 2940 : }
1005 :
1006 2940 : }
1007 :
1008 0 : void StokesImageUtil::directCFromR(ImageInterface<Complex>& out, const ImageInterface<Float>& in) {
1009 0 : AlwaysAssert(in.shape()(0)==out.shape()(0), AipsError);
1010 0 : AlwaysAssert(in.shape()(1)==out.shape()(1), AipsError);
1011 0 : AlwaysAssert(in.shape()(3)==out.shape()(3), AipsError);
1012 :
1013 0 : Vector<Int> inmap;
1014 0 : AlwaysAssert(StokesMap(inmap, in.coordinates()), AipsError);
1015 0 : Vector<Int> outmap;
1016 0 : AlwaysAssert(StokesMap(outmap, out.coordinates()), AipsError);
1017 :
1018 0 : Int innpol = in.shape()(inmap(2));
1019 0 : Int outnpol = out.shape()(outmap(2));
1020 0 : if(innpol != outnpol){
1021 0 : throw(AipsError("Cannot convert directly between images of different polarization shape"));
1022 : }
1023 0 : Vector<Int> inMap(innpol), outMap(outnpol);
1024 : StokesImageUtil::PolRep outPolFrame;
1025 0 : Int nStokesOut=CStokesPolMap(outMap, outPolFrame, out.coordinates());
1026 : StokesImageUtil::PolRep inPolFrame;
1027 0 : Int nStokesIn=CStokesPolMap(inMap, inPolFrame, in.coordinates());
1028 0 : AlwaysAssert(nStokesOut, AipsError);
1029 0 : AlwaysAssert(nStokesIn, AipsError);
1030 :
1031 0 : if(inPolFrame != outPolFrame){
1032 0 : throw(AipsError("Cannot convert directly between polarization types"));
1033 : }
1034 0 : IPosition inblc(4,0,0,0,0);
1035 0 : IPosition intrc(4,0,0,0,0);
1036 0 : intrc(inmap(0))=in.shape()(inmap(0))-1;
1037 0 : intrc(inmap(1))=in.shape()(inmap(1))-1;
1038 0 : intrc(inmap(2))=0;
1039 0 : intrc(inmap(3))= in.shape()(inmap(3))-1;
1040 0 : IPosition outblc(4,0,0,0,0);
1041 0 : IPosition outtrc(4,0,0,0,0);
1042 0 : outtrc(outmap(0))=in.shape()(outmap(0))-1;
1043 0 : outtrc(outmap(1))=in.shape()(outmap(1))-1;
1044 0 : outtrc(outmap(2))=0;
1045 0 : outtrc(outmap(3))= in.shape()(outmap(3))-1;
1046 :
1047 0 : for (Int k=0; k < innpol ; ++k){
1048 0 : inblc(inmap(2))=k;
1049 0 : intrc(inmap(2))=k;
1050 0 : Int outindex=-1;
1051 0 : for ( Int j=0; j < innpol; ++j){
1052 0 : if(inMap[k]==outMap[j])
1053 0 : outindex=j;
1054 : }
1055 0 : if(outindex < 0){
1056 0 : throw(AipsError("cannot match polarization in direct copy"));
1057 : }
1058 0 : outblc(outmap(2))=outindex;
1059 0 : outtrc(outmap(2))=outindex;
1060 0 : Slicer slin(inblc, intrc, Slicer::endIsLast);
1061 0 : Slicer slout(outblc, outtrc, Slicer::endIsLast);
1062 0 : SubImage<Complex> sliceout(out, slout, true);
1063 0 : SubImage<Float> slicein(in, slin);
1064 0 : sliceout.copyData(LatticeExpr<Complex>(toComplex(slicein)));
1065 0 : }
1066 :
1067 :
1068 :
1069 :
1070 0 : }
1071 :
1072 12 : void StokesImageUtil::directCToR(ImageInterface<Float>& out, ImageInterface<Complex>& in) {
1073 12 : AlwaysAssert(in.shape()(0)==out.shape()(0), AipsError);
1074 12 : AlwaysAssert(in.shape()(1)==out.shape()(1), AipsError);
1075 12 : AlwaysAssert(in.shape()(3)==out.shape()(3), AipsError);
1076 :
1077 12 : Vector<Int> inmap;
1078 12 : AlwaysAssert(StokesMap(inmap, in.coordinates()), AipsError);
1079 12 : Vector<Int> outmap;
1080 12 : AlwaysAssert(StokesMap(outmap, out.coordinates()), AipsError);
1081 :
1082 12 : Int innpol = in.shape()(inmap(2));
1083 12 : Int outnpol = out.shape()(outmap(2));
1084 12 : if(innpol != outnpol){
1085 0 : throw(AipsError("Cannot convert directly between images of different polarization shape"));
1086 : }
1087 12 : Vector<Int> inMap(innpol), outMap(outnpol);
1088 : StokesImageUtil::PolRep outPolFrame;
1089 12 : Int nStokesOut=CStokesPolMap(outMap, outPolFrame, out.coordinates());
1090 : StokesImageUtil::PolRep inPolFrame;
1091 12 : Int nStokesIn=CStokesPolMap(inMap, inPolFrame, in.coordinates());
1092 :
1093 :
1094 12 : AlwaysAssert(nStokesOut, AipsError);
1095 12 : AlwaysAssert(nStokesIn, AipsError);
1096 12 : if(inPolFrame != outPolFrame){
1097 0 : throw(AipsError("Cannot convert directly between polarization types"));
1098 : }
1099 12 : IPosition inblc(4,0,0,0,0);
1100 12 : IPosition intrc(4,0,0,0,0);
1101 12 : intrc(inmap(0))=in.shape()(inmap(0))-1;
1102 12 : intrc(inmap(1))=in.shape()(inmap(1))-1;
1103 12 : intrc(inmap(2))=0;
1104 12 : intrc(inmap(3))= in.shape()(inmap(3))-1;
1105 12 : IPosition outblc(4,0,0,0,0);
1106 12 : IPosition outtrc(4,0,0,0,0);
1107 12 : outtrc(outmap(0))=in.shape()(outmap(0))-1;
1108 12 : outtrc(outmap(1))=in.shape()(outmap(1))-1;
1109 12 : outtrc(outmap(2))=0;
1110 12 : outtrc(outmap(3))= in.shape()(outmap(3))-1;
1111 :
1112 30 : for (Int k=0; k < innpol ; ++k){
1113 18 : inblc(inmap(2))=k;
1114 18 : intrc(inmap(2))=k;
1115 18 : Int outindex=-1;
1116 48 : for ( Int j=0; j < innpol; ++j){
1117 30 : if(inMap[k]==outMap[j])
1118 18 : outindex=j;
1119 : }
1120 18 : if(outindex < 0){
1121 0 : throw(AipsError("cannot match polarization in direct copy"));
1122 : }
1123 18 : outblc(outmap(2))=outindex;
1124 18 : outtrc(outmap(2))=outindex;
1125 18 : Slicer slin(inblc, intrc, Slicer::endIsLast);
1126 18 : Slicer slout(outblc, outtrc, Slicer::endIsLast);
1127 18 : SubImage<Float> sliceout(out, slout, true);
1128 18 : SubImage<Complex> slicein(in, slin);
1129 18 : sliceout.copyData(LatticeExpr<Float>(real(slicein)));
1130 18 : }
1131 :
1132 :
1133 12 : }
1134 :
1135 2869 : void StokesImageUtil::To(ImageInterface<Float>& out, ImageInterface<Complex>& in) {
1136 :
1137 2869 : AlwaysAssert(in.shape()(0)==out.shape()(0), AipsError);
1138 2869 : AlwaysAssert(in.shape()(1)==out.shape()(1), AipsError);
1139 2869 : AlwaysAssert(in.shape()(3)==out.shape()(3), AipsError);
1140 :
1141 2869 : Vector<Int> map;
1142 2869 : AlwaysAssert(StokesMap(map, in.coordinates()), AipsError);
1143 2869 : Vector<Int> outmap;
1144 2869 : AlwaysAssert(StokesMap(outmap, out.coordinates()), AipsError);
1145 :
1146 2869 : Int nx = in.shape()(map(0));
1147 2869 : Int innpol = in.shape()(map(2));
1148 2869 : Int outnpol = out.shape()(outmap(2));
1149 :
1150 : // Loop over all planes
1151 5738 : LatticeStepper inls(in.shape(), IPosition(4, nx, 1, innpol, 1),
1152 5738 : IPosition(4, map(0), map(2), map(1), map(3)));
1153 5738 : LatticeStepper outls(out.shape(), IPosition(4, nx, 1, outnpol, 1),
1154 5738 : IPosition(4, map(0), map(2), map(1), map(3)));
1155 2869 : RO_LatticeIterator<Complex> inli(in, inls);
1156 2869 : LatticeIterator<Float> outli(out, outls);
1157 :
1158 2869 : Vector<Int> inMap(innpol), outMap(outnpol);
1159 2869 : StokesImageUtil::PolRep polFrame=StokesImageUtil::CIRCULAR;
1160 2869 : Int nStokesIn=CStokesPolMap(inMap, polFrame, in.coordinates());
1161 2869 : Int nStokesOut=StokesPolMap(outMap, out.coordinates());
1162 :
1163 2869 : if(nStokesOut <=0){
1164 12 : directCToR(out, in);
1165 12 : return;
1166 : }
1167 2857 : AlwaysAssert(nStokesOut, AipsError);
1168 : // Try taking real part only (uses LatticeExpr)
1169 2857 : if(nStokesIn==0) {
1170 2730 : nStokesIn=StokesPolMap(inMap, in.coordinates());
1171 2730 : if(nStokesIn==nStokesOut) {
1172 2730 : outli=LatticeIterator<Float> (out, inls);
1173 : //replaced copyData with iteration as it seems to load the whole array in memory
1174 3081664 : for (inli.reset(), outli.reset(); !inli.atEnd() && !outli.atEnd();
1175 3078934 : inli++,outli++) {
1176 3078934 : outli.woCursor()=real(inli.cursor());
1177 :
1178 : }
1179 :
1180 :
1181 : // out.copyData(LatticeExpr<Float>(real(in)));
1182 2730 : return;
1183 : }
1184 0 : throw(AipsError("Illegal conversion in ToStokesImage"));
1185 : }
1186 :
1187 127 : Vector<Complex> cs(4);
1188 127 : cs=Complex(0.0);
1189 127 : CStokesVector csv(0.0, 0.0, 0.0, 0.0);
1190 127 : StokesVector sv(0.0, 0.0, 0.0, 0.0);
1191 : Int pol;
1192 127 : if(nStokesIn==1) {
1193 0 : for (inli.reset(), outli.reset(); !inli.atEnd() && !outli.atEnd();
1194 0 : inli++,outli++) {
1195 0 : for (Int ix=0;ix<nx;ix++) {
1196 0 : cs(inMap(0))=inli.vectorCursor()(ix);
1197 0 : csv=CStokesVector(cs(0), cs(1), cs(2), cs(3));
1198 0 : if(polFrame==StokesImageUtil::LINEAR) {
1199 0 : applySlinInv(sv, csv);
1200 : }
1201 : else {
1202 0 : applyScircInv(sv, csv);
1203 : }
1204 0 : if(nStokesOut==1) {
1205 0 : outli.rwVectorCursor()(ix)=sv(outMap(0));
1206 : }
1207 : else {
1208 0 : for (pol=0;pol<outnpol;pol++)
1209 0 : outli.rwMatrixCursor()(ix,pol)=sv(outMap(pol));
1210 : }
1211 : }
1212 : }
1213 : }
1214 : else {
1215 24027 : for (inli.reset(), outli.reset(); !inli.atEnd() && !outli.atEnd(); inli++,outli++) {
1216 2413900 : for (Int ix=0;ix<nx;ix++) {
1217 2390000 : cs=Complex(0.0);
1218 9970000 : for (pol=0;pol<innpol;pol++) cs(inMap(pol))=inli.matrixCursor()(ix,pol);
1219 2390000 : csv=CStokesVector(cs(0), cs(1), cs(2), cs(3));
1220 2390000 : if(polFrame==StokesImageUtil::LINEAR) {
1221 860000 : applySlinInv(sv, csv);
1222 : }
1223 : else {
1224 1530000 : applyScircInv(sv, csv);
1225 : }
1226 2390000 : if(nStokesOut==1) {
1227 830000 : outli.rwVectorCursor()(ix)=sv(outMap(0));
1228 : }
1229 : else {
1230 7480000 : for (pol=0;pol<outnpol;pol++) {
1231 5920000 : outli.rwMatrixCursor()(ix,pol)=sv(outMap(pol));
1232 : }
1233 : }
1234 : }
1235 : }
1236 : }
1237 22063 : };
1238 :
1239 1163 : void StokesImageUtil::ToStokesPSF(ImageInterface<Float>& out, ImageInterface<Complex>& in) {
1240 :
1241 1163 : AlwaysAssert(in.shape()(0)==out.shape()(0), AipsError);
1242 1163 : AlwaysAssert(in.shape()(1)==out.shape()(1), AipsError);
1243 1163 : AlwaysAssert(in.shape()(3)==out.shape()(3), AipsError);
1244 :
1245 1163 : Vector<Int> map;
1246 1163 : AlwaysAssert(StokesMap(map, in.coordinates()), AipsError);
1247 1163 : Vector<Int> outmap;
1248 1163 : AlwaysAssert(StokesMap(outmap, out.coordinates()), AipsError);
1249 :
1250 1163 : Int nx = in.shape()(map(0));
1251 1163 : Int innpol = in.shape()(map(2));
1252 1163 : Int outnpol = out.shape()(outmap(2));
1253 :
1254 : // Loop over all planes
1255 2326 : LatticeStepper inls(in.shape(), IPosition(4, nx, 1, innpol, 1),
1256 2326 : IPosition(4, map(0), map(2), map(1), map(3)));
1257 2326 : LatticeStepper outls(out.shape(), IPosition(4, nx, 1, outnpol, 1),
1258 2326 : IPosition(4, map(0), map(2), map(1), map(3)));
1259 1163 : RO_LatticeIterator<Complex> inli(in, inls);
1260 1163 : LatticeIterator<Float> outli(out, outls);
1261 :
1262 1163 : Vector<Int> inMap(innpol), outMap(outnpol);
1263 1163 : StokesImageUtil::PolRep polFrame=StokesImageUtil::CIRCULAR;
1264 1163 : Int nStokesIn=CStokesPolMap(inMap, polFrame, in.coordinates());
1265 1163 : Int nStokesOut=StokesPolMap(outMap, out.coordinates());
1266 1163 : if(nStokesOut <=0){
1267 0 : directCToR(out, in);
1268 0 : return;
1269 : }
1270 1163 : AlwaysAssert(nStokesOut, AipsError);
1271 : // Try taking real part only (uses LatticeExpr)
1272 1163 : if(nStokesIn==0) {
1273 1093 : nStokesIn=StokesPolMap(inMap, in.coordinates());
1274 1093 : if(nStokesIn==nStokesOut) {
1275 : //replaced copyData with iteration as it seems to load the whole array in memory
1276 1093 : outli=LatticeIterator<Float> (out, inls);
1277 1184389 : for (inli.reset(), outli.reset(); !inli.atEnd() && !outli.atEnd();
1278 1183296 : inli++,outli++) {
1279 1183296 : outli.woCursor()=real(inli.cursor());
1280 :
1281 : }
1282 :
1283 : //out.copyData(LatticeExpr<Float>(real(in)));
1284 1093 : return;
1285 : }
1286 0 : throw(AipsError("Illegal conversion in ToStokesImage"));
1287 : }
1288 :
1289 : /* STOKESDBG */ //cout << "ToStokesPSF - nstokesIn : " << nStokesIn << " inmap : " << inMap << " nstokesOut : " << nStokesOut << " outmap : " << outMap << endl;
1290 : /* STOKESDBG */ //cout << "IN (data) stokesCoord : " << ( (in.coordinates()).stokesCoordinate( (in.coordinates()).findCoordinate(Coordinate::STOKES) ) ).stokes() << " OUT (image) stokesCoord : " << ( (out.coordinates()).stokesCoordinate( (out.coordinates()).findCoordinate(Coordinate::STOKES) ) ).stokes() << endl;
1291 :
1292 70 : Vector<Int> inST = ( (in.coordinates()).stokesCoordinate( (in.coordinates()).findCoordinate(Coordinate::STOKES) ) ).stokes();
1293 70 : Vector<Int> outST = ( (out.coordinates()).stokesCoordinate( (out.coordinates()).findCoordinate(Coordinate::STOKES) ) ).stokes();
1294 :
1295 70 : Vector<Complex> cs(4);
1296 70 : cs=Complex(0.0);
1297 70 : CStokesVector csv(0.0, 0.0, 0.0, 0.0);
1298 70 : StokesVector sv(0.0, 0.0, 0.0, 0.0);
1299 : Int pol;
1300 70 : if(nStokesIn==1) {
1301 0 : for (inli.reset(), outli.reset(); !inli.atEnd() && !outli.atEnd();
1302 0 : inli++,outli++) {
1303 0 : for (Int ix=0;ix<nx;ix++) {
1304 0 : cs(inMap(0))=inli.vectorCursor()(ix);
1305 0 : csv=CStokesVector(cs(0), cs(1), cs(2), cs(3));
1306 0 : if(polFrame==StokesImageUtil::LINEAR) {
1307 0 : applySlinInv(sv, csv);
1308 : }
1309 : else {
1310 0 : applyScircInv(sv, csv);
1311 : }
1312 : // nstokesout will be 1, 2 or 4 only.
1313 0 : if(nStokesOut==1) {
1314 : //outli.rwVectorCursor()(ix)=sv(0); //sv(outMap(0));
1315 : // If outstokescoord = I or Q, use (I:0) : XX+YY for the PSF
1316 : // If outstokescoord = U or V, use (U:2) : XY+YX for the PSF
1317 0 : if( polFrame==StokesImageUtil::LINEAR )
1318 : {
1319 0 : if(Stokes::type(outST[0])==Stokes::I || Stokes::type(outST[0])==Stokes::Q )
1320 0 : outli.rwVectorCursor()(ix)=sv(0);
1321 0 : if( Stokes::type(outST[0])==Stokes::U || Stokes::type(outST[0])==Stokes::V )
1322 0 : outli.rwVectorCursor()(ix)=sv(2);
1323 : }
1324 0 : if( polFrame==StokesImageUtil::CIRCULAR )
1325 : {
1326 0 : if(Stokes::type(outST[0])==Stokes::I || Stokes::type(outST[0])==Stokes::V )
1327 0 : outli.rwVectorCursor()(ix)=sv(0);
1328 0 : if( Stokes::type(outST[0])==Stokes::Q || Stokes::type(outST[0])==Stokes::U )
1329 0 : outli.rwVectorCursor()(ix)=sv(1);
1330 : }
1331 : }
1332 0 : else if(nStokesOut==2) {
1333 : // If outstokescoord = IQ, use (I:0) : XX+YY for the PSF in both planes
1334 : // If outstokescoord = UV, use (U:2) : XY+YX for the PSF in both planes
1335 0 : if( polFrame==StokesImageUtil::LINEAR )
1336 : {
1337 0 : if(Stokes::type(outST[0])==Stokes::I || Stokes::type(outST[0])==Stokes::Q ||
1338 0 : Stokes::type(outST[1])==Stokes::I || Stokes::type(outST[1])==Stokes::Q )
1339 0 : { outli.rwMatrixCursor()(ix,0)=sv(0); outli.rwMatrixCursor()(ix,1)=sv(0); }
1340 0 : if( Stokes::type(outST[0])==Stokes::U || Stokes::type(outST[0])==Stokes::V ||
1341 0 : Stokes::type(outST[1])==Stokes::U || Stokes::type(outST[1])==Stokes::V )
1342 0 : { outli.rwMatrixCursor()(ix,0)=sv(2); outli.rwMatrixCursor()(ix,1)=sv(2); }
1343 : }
1344 0 : if( polFrame==StokesImageUtil::CIRCULAR )
1345 : {
1346 0 : if(Stokes::type(outST[0])==Stokes::I || Stokes::type(outST[0])==Stokes::V ||
1347 0 : Stokes::type(outST[1])==Stokes::I || Stokes::type(outST[1])==Stokes::V )
1348 0 : { outli.rwMatrixCursor()(ix,0)=sv(0); outli.rwMatrixCursor()(ix,1)=sv(0); }
1349 0 : if( Stokes::type(outST[0])==Stokes::Q || Stokes::type(outST[0])==Stokes::U ||
1350 0 : Stokes::type(outST[1])==Stokes::Q || Stokes::type(outST[1])==Stokes::U )
1351 0 : { outli.rwMatrixCursor()(ix,0)=sv(1); outli.rwMatrixCursor()(ix,1)=sv(1); }
1352 : }
1353 : }
1354 : else { // nstokesout = 4 : do a one-to-one mapping.
1355 0 : for (pol=0;pol<outnpol;pol++)
1356 0 : outli.rwMatrixCursor()(ix,pol)=sv(outMap(pol));
1357 : }
1358 : }
1359 : }
1360 : }
1361 : else {
1362 12770 : for (inli.reset(), outli.reset(); !inli.atEnd() && !outli.atEnd(); inli++,outli++) {
1363 1282700 : for (Int ix=0;ix<nx;ix++) {
1364 1270000 : cs=Complex(0.0);
1365 5270000 : for (pol=0;pol<innpol;pol++) cs(inMap(pol))=inli.matrixCursor()(ix,pol);
1366 1270000 : csv=CStokesVector(cs(0), cs(1), cs(2), cs(3));
1367 1270000 : if(polFrame==StokesImageUtil::LINEAR) {
1368 430000 : applySlinInv(sv, csv);
1369 : }
1370 : else {
1371 840000 : applyScircInv(sv, csv);
1372 : }
1373 1270000 : if(nStokesOut==1) {
1374 : //outli.rwVectorCursor()(ix)=sv(0); //sv(outMap(0));
1375 : // If outstokescoord = I or Q, use (I:0)
1376 : // If outstokescoord = U or V, use (U:2)
1377 460000 : if( polFrame==StokesImageUtil::LINEAR )
1378 : {
1379 0 : if(Stokes::type(outST[0])==Stokes::I || Stokes::type(outST[0])==Stokes::Q )
1380 0 : outli.rwVectorCursor()(ix)=sv(0);
1381 0 : if( Stokes::type(outST[0])==Stokes::U || Stokes::type(outST[0])==Stokes::V )
1382 0 : outli.rwVectorCursor()(ix)=sv(2);
1383 : }
1384 460000 : if( polFrame==StokesImageUtil::CIRCULAR )
1385 : {
1386 460000 : if(Stokes::type(outST[0])==Stokes::I || Stokes::type(outST[0])==Stokes::V )
1387 40000 : outli.rwVectorCursor()(ix)=sv(0);
1388 460000 : if( Stokes::type(outST[0])==Stokes::Q || Stokes::type(outST[0])==Stokes::U )
1389 420000 : outli.rwVectorCursor()(ix)=sv(1);
1390 : }
1391 : }
1392 810000 : else if(nStokesOut==2) {
1393 : // If outstokescoord = IQ, use (I:0) : XX+YY for the PSF in both planes
1394 : // If outstokescoord = UV, use (U:2) : XY+YX for the PSF in both planes
1395 80000 : if( polFrame==StokesImageUtil::LINEAR )
1396 : {
1397 0 : if(Stokes::type(outST[0])==Stokes::I || Stokes::type(outST[0])==Stokes::Q ||
1398 0 : Stokes::type(outST[1])==Stokes::I || Stokes::type(outST[1])==Stokes::Q )
1399 0 : { outli.rwMatrixCursor()(ix,0)=sv(0); outli.rwMatrixCursor()(ix,1)=sv(0); }
1400 0 : if( Stokes::type(outST[0])==Stokes::U || Stokes::type(outST[0])==Stokes::V ||
1401 0 : Stokes::type(outST[1])==Stokes::U || Stokes::type(outST[1])==Stokes::V )
1402 0 : { outli.rwMatrixCursor()(ix,0)=sv(2); outli.rwMatrixCursor()(ix,1)=sv(2); }
1403 : }
1404 80000 : if( polFrame==StokesImageUtil::CIRCULAR )
1405 : {
1406 120000 : if(Stokes::type(outST[0])==Stokes::I || Stokes::type(outST[0])==Stokes::V ||
1407 120000 : Stokes::type(outST[1])==Stokes::I || Stokes::type(outST[1])==Stokes::V )
1408 40000 : { outli.rwMatrixCursor()(ix,0)=sv(0); outli.rwMatrixCursor()(ix,1)=sv(0); }
1409 120000 : if( Stokes::type(outST[0])==Stokes::Q || Stokes::type(outST[0])==Stokes::U ||
1410 120000 : Stokes::type(outST[1])==Stokes::Q || Stokes::type(outST[1])==Stokes::U )
1411 40000 : { outli.rwMatrixCursor()(ix,0)=sv(1); outli.rwMatrixCursor()(ix,1)=sv(1); }
1412 : }
1413 : }
1414 : else { // nstokesout = 4
1415 3650000 : for (pol=0;pol<outnpol;pol++) {
1416 2920000 : outli.rwMatrixCursor()(ix,pol)=sv(outMap(0));
1417 : }
1418 : }
1419 : }
1420 : }
1421 : }
1422 8814 : };
1423 :
1424 747 : void StokesImageUtil::ToStokesSumWt(Matrix<Float> sumwtStokes, Matrix<Float> sumwtCorr, const CoordinateSystem& outcoord, const CoordinateSystem& incoord)
1425 : {
1426 :
1427 747 : AlwaysAssert( sumwtStokes.shape()[1] == sumwtCorr.shape()(1) , AipsError );
1428 747 : Vector<Int> inmap;
1429 747 : AlwaysAssert(StokesMap(inmap, incoord), AipsError);
1430 747 : Vector<Int> outmap;
1431 747 : AlwaysAssert(StokesMap(outmap, outcoord), AipsError);
1432 :
1433 :
1434 747 : Int innpol = sumwtCorr.shape()[0];
1435 747 : Int outnpol = sumwtStokes.shape()[0];
1436 :
1437 :
1438 747 : Vector<Int> inMap(innpol), outMap(outnpol);
1439 747 : StokesImageUtil::PolRep polFrame=StokesImageUtil::CIRCULAR;
1440 747 : Int nStokesIn=CStokesPolMap(inMap, polFrame, incoord);
1441 747 : Int nStokesOut=StokesPolMap(outMap, outcoord);
1442 747 : if(nStokesOut <=0 || nStokesIn==0) {
1443 684 : if(innpol != outnpol)
1444 0 : throw(AipsError("Cannot copy sumwt as input and output have different pol axis length"));
1445 684 : convertArray(sumwtStokes, sumwtCorr);
1446 684 : return;
1447 : }
1448 63 : if( nStokesOut==1 && nStokesIn==2){
1449 :
1450 79 : for(uInt chan=0;chan<sumwtStokes.shape()[1];chan++)
1451 : {
1452 50 : sumwtStokes(0,chan) = min(sumwtCorr(0,chan),sumwtCorr(1, chan));
1453 : }
1454 : }
1455 63 : if( nStokesOut==2 && nStokesIn==2){
1456 :
1457 12 : for(uInt chan=0;chan<sumwtStokes.shape()[1];chan++)
1458 : {
1459 8 : sumwtStokes(0,chan) = min(sumwtCorr(0,chan),sumwtCorr(1, chan));
1460 8 : sumwtStokes(1,chan) = min(sumwtCorr(1,chan),sumwtCorr(1, chan));
1461 : }
1462 : }
1463 63 : if( nStokesOut==4 && nStokesIn==4){
1464 :
1465 91 : for(uInt chan=0;chan<sumwtStokes.shape()[1];chan++)
1466 : {
1467 61 : if(polFrame==StokesImageUtil::LINEAR){
1468 37 : sumwtStokes(0,chan) = min(sumwtCorr(0,chan),sumwtCorr(3, chan));
1469 37 : sumwtStokes(1,chan) = min(sumwtCorr(0,chan),sumwtCorr(3, chan));
1470 37 : sumwtStokes(2,chan) = min(sumwtCorr(1,chan),sumwtCorr(2, chan));
1471 37 : sumwtStokes(3,chan) = min(sumwtCorr(1,chan),sumwtCorr(2, chan));
1472 : }
1473 : else{
1474 24 : sumwtStokes(0,chan) = min(sumwtCorr(0,chan),sumwtCorr(3, chan));
1475 24 : sumwtStokes(1,chan) = min(sumwtCorr(1,chan),sumwtCorr(2, chan));
1476 24 : sumwtStokes(2,chan) = min(sumwtCorr(1,chan),sumwtCorr(2, chan));
1477 24 : sumwtStokes(3,chan) = min(sumwtCorr(0,chan),sumwtCorr(3, chan));
1478 : }
1479 : }
1480 : }
1481 :
1482 :
1483 2799 : }
1484 :
1485 675 : void StokesImageUtil::ToStokesSumWt(Matrix<Float> sumwtStokes, Matrix<Float> sumwtCorr)
1486 : {
1487 675 : AlwaysAssert( sumwtStokes.shape()[1] == sumwtCorr.shape()(1) , AipsError ); //same nchan
1488 675 : AlwaysAssert( sumwtCorr.shape()[0] > 0, AipsError ); // at least one correlation gridded.
1489 :
1490 : /// Pick the value from the first correlation plane, and fill it into all stokes planes.
1491 : /// This is valid when the same weights are used across correlations.
1492 1386 : for(uInt pol=0;pol<sumwtStokes.shape()[0];pol++)
1493 : {
1494 1422 : for(uInt chan=0;chan<sumwtStokes.shape()[1];chan++)
1495 : {
1496 711 : sumwtStokes(pol,chan) = sumwtCorr(0,chan);
1497 : }
1498 : }
1499 :
1500 675 : }
1501 :
1502 :
1503 : #if 0
1504 : void StokesImageUtil::ToStokesPSF(ImageInterface<Float>& out, ImageInterface<Complex>& in) {
1505 :
1506 : AlwaysAssert(in.shape()(0)==out.shape()(0), AipsError);
1507 : AlwaysAssert(in.shape()(1)==out.shape()(1), AipsError);
1508 : AlwaysAssert(out.shape()(2)==out.shape()(3), AipsError);
1509 : AlwaysAssert(in.shape()(3)==out.shape()(4), AipsError);
1510 :
1511 : Vector<Int> map;
1512 : AlwaysAssert(StokesMap(map, in.coordinates()), AipsError);
1513 : Vector<Int> outmap;
1514 : AlwaysAssert(StokesMap(outmap, out.coordinates()), AipsError);
1515 :
1516 : Int nx = in.shape()(map(0));
1517 : Int innpol = in.shape()(map(2));
1518 : Int outnpol = out.shape()(outmap(2));
1519 :
1520 : // Loop over all planes. For the input we get a line of complex pixels
1521 : // for the output we put a line of float matrices
1522 : LatticeStepper inls(in.shape(), IPosition(4, nx, 1, innpol, 1),
1523 : IPosition(4, map(0), map(2), map(1), map(3)));
1524 : LatticeStepper outls(out.shape(), IPosition(5, nx, 1, outnpol, outnpol, 1),
1525 : IPosition(5, map(0), map(3), map(1), map(2), 4));
1526 : RO_LatticeIterator<Complex> inli(in, inls);
1527 : LatticeIterator<Float> outli(out, outls);
1528 :
1529 : Vector<Int> inMap(innpol), outMap(outnpol);
1530 : StokesImageUtil::PolRep polFrame;
1531 : Int nStokesIn=CStokesPolMap(inMap, polFrame, in.coordinates());
1532 : AlwaysAssert(nStokesIn, AipsError);
1533 : Int nStokesOut=StokesPolMap(outMap, out.coordinates());
1534 : AlwaysAssert(nStokesOut, AipsError);
1535 :
1536 : Matrix<Complex> s(4,4), sAdjoint(4,4);
1537 : s=0;
1538 : if(polFrame==StokesImageUtil::LINEAR) {
1539 : s(0,0)=Complex(0.5);
1540 : s(0,1)=Complex(0.5);
1541 : s(1,2)=Complex(0.5);
1542 : s(1,3)=Complex(0.0,0.5);
1543 : s(2,2)=Complex(0.5);
1544 : s(2,3)=Complex(0.0,-0.5);
1545 : s(3,0)=Complex(0.5);
1546 : s(3,1)=Complex(-0.5);
1547 : }
1548 : else {
1549 : s(0,0)=Complex(0.5);
1550 : s(0,3)=Complex(0.5);
1551 : s(1,1)=Complex(0.5);
1552 : s(1,2)=Complex(0.0,0.5);
1553 : s(2,1)=Complex(0.5);
1554 : s(2,2)=Complex(0.0,-0.5);
1555 : s(3,0)=Complex(0.5);
1556 : s(3,3)=Complex(-0.5);
1557 : }
1558 : sAdjoint=adjoint(s);
1559 :
1560 : Matrix<Complex> lambda(4,4);
1561 : lambda=Complex(0.0);
1562 :
1563 : Matrix<Float> stokesPSF(4,4);
1564 : Int i, j;
1565 : for (inli.reset(), outli.reset(); !inli.atEnd() && !outli.atEnd(); inli++,outli++) {
1566 : for (Int ix=0;ix<nx;ix++) {
1567 : // Get the channel
1568 : // Fill in the weight vector
1569 : if(innpol==1) {
1570 : lambda(inMap(0),inMap(0))=inli.vectorCursor()(ix);
1571 : }
1572 : else {
1573 : for (i=0;i<innpol;i++) lambda(inMap(i),inMap(i))=
1574 : inli.matrixCursor()(ix,i);
1575 : }
1576 : stokesPSF=real(product(s,product(lambda, sAdjoint)));
1577 : if(nStokesOut==1) {
1578 : outli.rwVectorCursor()(ix)=stokesPSF(outMap(0),outMap(0));
1579 : }
1580 : else {
1581 : for (j=0;j<outnpol;j++) {
1582 : for (i=0;i<outnpol;i++) {
1583 : outli.rwCubeCursor()(ix,i,j)=stokesPSF(outMap(i),outMap(j));
1584 : }
1585 : }
1586 : }
1587 : }
1588 : }
1589 : };
1590 : #endif
1591 :
1592 1307 : void StokesImageUtil::From(ImageInterface<Complex>& out,
1593 : const ImageInterface<Float>& in) {
1594 :
1595 1307 : AlwaysAssert(in.shape()(0)==out.shape()(0), AipsError);
1596 1307 : AlwaysAssert(in.shape()(1)==out.shape()(1), AipsError);
1597 1307 : AlwaysAssert(in.shape()(3)==out.shape()(3), AipsError);
1598 :
1599 :
1600 1307 : Vector<Int> map;
1601 1307 : AlwaysAssert(StokesMap(map, in.coordinates()), AipsError);
1602 1307 : Vector<Int> outmap;
1603 1307 : AlwaysAssert(StokesMap(outmap, out.coordinates()), AipsError);
1604 :
1605 1307 : Int nx = in.shape()(map(0));
1606 1307 : Int innpol = in.shape()(map(2));
1607 1307 : Int outnpol = out.shape()(outmap(2));
1608 :
1609 1307 : Vector<Int> inMap(innpol), outMap(outnpol);
1610 1307 : StokesImageUtil::PolRep polFrame=StokesImageUtil::LINEAR;
1611 1307 : Int nStokesIn=StokesPolMap(inMap, in.coordinates());
1612 1307 : if(nStokesIn <=0){
1613 0 : directCFromR(out, in);
1614 0 : return;
1615 : }
1616 1307 : AlwaysAssert(nStokesIn, AipsError);
1617 1307 : Int nStokesOut=CStokesPolMap(outMap, polFrame, out.coordinates());
1618 1307 : if(nStokesOut==0) {
1619 1244 : nStokesOut=StokesPolMap(outMap, out.coordinates());
1620 1244 : if(nStokesIn==nStokesOut) {
1621 1244 : out.copyData(LatticeExpr<Complex>(toComplex(in)));
1622 1244 : return;
1623 : }
1624 0 : throw(AipsError("Illegal conversion in FromStokesImage"));
1625 : }
1626 :
1627 : // Loop over all planes
1628 126 : LatticeStepper inls(in.shape(), IPosition(4, nx, 1, innpol, 1),
1629 126 : IPosition(4, map(0), map(2), map(1), map(3)));
1630 126 : LatticeStepper outls(out.shape(), IPosition(4, nx, 1, outnpol, 1),
1631 126 : IPosition(4, map(0), map(2), map(1), map(3)));
1632 63 : RO_LatticeIterator<Float> inli(in, inls);
1633 63 : LatticeIterator<Complex> outli(out, outls);
1634 :
1635 63 : Vector<Float> s(4);
1636 63 : CStokesVector csv(0.0, 0.0, 0.0, 0.0);
1637 63 : StokesVector sv(0.0, 0.0, 0.0, 0.0);
1638 : Int pol;
1639 63 : if(nStokesIn==1) {
1640 3821 : for (inli.reset(), outli.reset(); !inli.atEnd() && !outli.atEnd(); inli++,outli++) {
1641 383800 : for (Int ix=0;ix<nx;ix++) {
1642 380000 : s=0.0;
1643 380000 : s(inMap(0))=inli.vectorCursor()(ix);
1644 380000 : sv=StokesVector(s(0), s(1), s(2), s(3));
1645 380000 : if(polFrame==StokesImageUtil::LINEAR) {
1646 0 : applySlin(csv, sv);
1647 : }
1648 : else {
1649 380000 : applyScirc(csv, sv);
1650 : }
1651 380000 : if(nStokesOut==1) {
1652 0 : outli.rwVectorCursor()(ix)=csv(outMap(0));
1653 : }
1654 : else {
1655 1140000 : for (pol=0;pol<outnpol;pol++) {
1656 760000 : outli.rwMatrixCursor()(ix,pol)=csv(outMap(pol));
1657 : }
1658 : }
1659 : }
1660 : }
1661 : }
1662 : else {
1663 8342 : for (inli.reset(), outli.reset(); !inli.atEnd() && !outli.atEnd(); inli++,outli++) {
1664 838300 : for (Int ix=0;ix<nx;ix++) {
1665 830000 : s=0.0;
1666 3990000 : for(pol=0;pol<innpol;pol++) s(inMap(pol))=inli.matrixCursor()(ix,pol);
1667 830000 : sv=StokesVector(s(0), s(1), s(2), s(3));
1668 830000 : if(polFrame==StokesImageUtil::LINEAR) {
1669 490000 : applySlin(csv, sv);
1670 : }
1671 : else {
1672 340000 : applyScirc(csv, sv);
1673 : }
1674 830000 : if(nStokesOut==1) {
1675 0 : outli.rwVectorCursor()(ix)=csv(outMap(0));
1676 : }
1677 : else {
1678 3990000 : for (pol=0;pol<outnpol;pol++) {
1679 3160000 : outli.rwMatrixCursor()(ix,pol)=csv(outMap(pol));
1680 : }
1681 : }
1682 : }
1683 : }
1684 : }
1685 5039 : };
1686 :
1687 : // map(axis) is the polarization in the sequence I,Q,U,V
1688 : // This is the sequence used by StokesVector
1689 11153 : Int StokesImageUtil::StokesPolMap(Vector<Int>& map, const CoordinateSystem& coord) {
1690 :
1691 11153 : map=-1;
1692 : Int stokesIndex;
1693 11153 : stokesIndex= coord.findCoordinate(Coordinate::STOKES);
1694 11153 : StokesCoordinate stokesCoord=coord.stokesCoordinate(stokesIndex);
1695 : Int pol, p;
1696 11153 : Bool Found= stokesCoord.toPixel(p, Stokes::I)||
1697 212 : stokesCoord.toPixel(p, Stokes::Q)||
1698 11395 : stokesCoord.toPixel(p, Stokes::U)||
1699 30 : stokesCoord.toPixel(p, Stokes::V);
1700 11153 : if(Found) {
1701 11137 : pol=0;
1702 11137 : Int found=0;
1703 11137 : if(stokesCoord.toPixel(p, Stokes::I)) {map(p)=pol;found++;} pol++;
1704 11137 : if(stokesCoord.toPixel(p, Stokes::Q)) {map(p)=pol;found++;} pol++;
1705 11137 : if(stokesCoord.toPixel(p, Stokes::U)) {map(p)=pol;found++;} pol++;
1706 11137 : if(stokesCoord.toPixel(p, Stokes::V)) {map(p)=pol;found++;} pol++;
1707 11137 : return found;
1708 : }
1709 16 : return 0;
1710 11153 : }
1711 :
1712 : // map(axis) is the polarization in the sequence XX,XY,YX,YY or RR,RL,LR,LL.
1713 : // This is the sequence used by CStokesVector
1714 6110 : Int StokesImageUtil::CStokesPolMap(Vector<Int>& map, StokesImageUtil::PolRep& polFrame,
1715 : const CoordinateSystem& coord) {
1716 :
1717 6110 : map=-1;
1718 : Int stokesIndex;
1719 6110 : stokesIndex=coord.findCoordinate(Coordinate::STOKES);
1720 6110 : StokesCoordinate stokesCoord=coord.stokesCoordinate(stokesIndex);
1721 : Int pol, p;
1722 6110 : Int found=0;
1723 6110 : Bool Linear= stokesCoord.toPixel(p, Stokes::XX)||
1724 5965 : stokesCoord.toPixel(p, Stokes::XY)||
1725 18040 : stokesCoord.toPixel(p, Stokes::YX)||
1726 5965 : stokesCoord.toPixel(p, Stokes::YY);
1727 6110 : if(Linear) {
1728 155 : pol=0;
1729 155 : if(stokesCoord.toPixel(p, Stokes::XX)) {map(p)=pol;found++;} pol++;
1730 155 : if(stokesCoord.toPixel(p, Stokes::XY)) {map(p)=pol;found++;} pol++;
1731 155 : if(stokesCoord.toPixel(p, Stokes::YX)) {map(p)=pol;found++;} pol++;
1732 155 : if(stokesCoord.toPixel(p, Stokes::YY)) {map(p)=pol;found++;} pol++;
1733 155 : polFrame=StokesImageUtil::LINEAR;
1734 155 : return found;
1735 : }
1736 5955 : Bool Circular= stokesCoord.toPixel(p, Stokes::LL)||
1737 5871 : stokesCoord.toPixel(p, Stokes::LR)||
1738 17573 : stokesCoord.toPixel(p, Stokes::RL)||
1739 5747 : stokesCoord.toPixel(p, Stokes::RR);
1740 :
1741 :
1742 5955 : if(Circular) {
1743 208 : pol=0;
1744 208 : if(stokesCoord.toPixel(p, Stokes::RR)) {map(p)=pol;found++;} pol++;
1745 208 : if(stokesCoord.toPixel(p, Stokes::RL)) {map(p)=pol;found++;} pol++;
1746 208 : if(stokesCoord.toPixel(p, Stokes::LR)) {map(p)=pol;found++;} pol++;
1747 208 : if(stokesCoord.toPixel(p, Stokes::LL)) {map(p)=pol;found++;} pol++;
1748 208 : polFrame=StokesImageUtil::CIRCULAR;
1749 208 : return found;
1750 : }
1751 5747 : return 0;
1752 6110 : }
1753 :
1754 0 : CoordinateSystem StokesImageUtil::StokesCoordFromMS(const IPosition& shape, Vector<Double>& deltas,
1755 : MeasurementSet& ms) {
1756 0 : Vector<Int> whichStokes(0);
1757 0 : return StokesCoordFromMS(shape, deltas, ms, whichStokes);
1758 0 : }
1759 :
1760 0 : CoordinateSystem StokesImageUtil::StokesCoordFromMS(const IPosition& shape, Vector<Double>& deltas,
1761 : MeasurementSet& ms, Vector<Int>& whichStokes,
1762 : Bool doCStokes, Int fieldID, Int SPWID, Int feedID) {
1763 :
1764 0 : MSColumns msc(ms);
1765 :
1766 0 : Int nx=shape(0);
1767 0 : Int ny=shape(1);
1768 0 : Int npol=shape(2);
1769 0 : Int nchan=shape(3);
1770 :
1771 0 : Vector<Double> refCoord=msc.field().phaseDir()(fieldID);
1772 0 : Vector<Double> refPixel(2);
1773 0 : refPixel(0)=Double(nx/2);
1774 0 : refPixel(1)=Double(ny/2);
1775 :
1776 0 : Matrix<Double> xform(2,2);
1777 0 : xform=0.0;xform.diagonal()=1.0;
1778 : DirectionCoordinate myRaDec(MDirection::J2000,
1779 : Projection::SIN,
1780 0 : refCoord(0), refCoord(1),
1781 0 : deltas(0), deltas(1),
1782 : xform,
1783 0 : refPixel(0), refPixel(1));
1784 :
1785 : // Frequency
1786 0 : Vector<Double> chanFreq=msc.spectralWindow().chanFreq()(SPWID);
1787 0 : if(nchan==0) nchan=chanFreq.shape()(SPWID);
1788 0 : Double refChan=0.0;
1789 0 : Vector<Double> freqResolution=msc.spectralWindow().resolution()(SPWID);
1790 :
1791 : // Retrieve the first rest frequency used with this SPW_ID (for now)
1792 0 : MSDopplerUtil msdoppler(ms);
1793 0 : Vector<Double> restFreqArray;
1794 : Double restFreq;
1795 0 : if (msdoppler.dopplerInfo(restFreqArray,SPWID,fieldID)) {
1796 0 : restFreq=restFreqArray(0);
1797 : } else {
1798 0 : restFreq=1.0;
1799 : };
1800 :
1801 0 : SpectralCoordinate mySpectral(MFrequency::TOPO, chanFreq(0),
1802 0 : freqResolution(0), refChan,
1803 0 : restFreq);
1804 :
1805 : // Polarization: If the specified whichStokes are ok, we use them
1806 : // otherwise we guess
1807 0 : if(Int(whichStokes.nelements())!=npol) {
1808 0 : Vector<String> polType=msc.feed().polarizationType()(feedID);
1809 0 : whichStokes.resize(npol);
1810 : // Polarization
1811 0 : if(doCStokes) {
1812 0 : if (polType(0)=="X" || polType(0)=="Y") {
1813 0 : switch(npol) {
1814 0 : case 1:
1815 0 : whichStokes.resize(1);
1816 0 : whichStokes(0)=Stokes::I;
1817 0 : break;
1818 0 : case 2:
1819 0 : whichStokes.resize(2);
1820 0 : whichStokes(0)=Stokes::XX;
1821 0 : whichStokes(1)=Stokes::YY;
1822 0 : break;
1823 0 : default:
1824 0 : whichStokes.resize(4);
1825 0 : whichStokes(0)=Stokes::XX;
1826 0 : whichStokes(1)=Stokes::XY;
1827 0 : whichStokes(2)=Stokes::YX;
1828 0 : whichStokes(3)=Stokes::YY;
1829 : }
1830 : }
1831 : else {
1832 0 : switch(npol) {
1833 0 : case 1:
1834 0 : whichStokes.resize(1);
1835 0 : whichStokes(0)=Stokes::I;
1836 0 : break;
1837 0 : case 2:
1838 0 : whichStokes.resize(2);
1839 0 : whichStokes(0)=Stokes::LL;
1840 0 : whichStokes(1)=Stokes::RR;
1841 0 : break;
1842 0 : default:
1843 0 : whichStokes.resize(4);
1844 0 : whichStokes(0)=Stokes::LL;
1845 0 : whichStokes(1)=Stokes::LR;
1846 0 : whichStokes(2)=Stokes::RL;
1847 0 : whichStokes(3)=Stokes::RR;
1848 : };
1849 : }
1850 : }
1851 : else {
1852 0 : switch(npol) {
1853 0 : case 1:
1854 0 : whichStokes.resize(1);
1855 0 : whichStokes(0)=Stokes::I;
1856 0 : break;
1857 0 : case 2:
1858 0 : whichStokes.resize(2);
1859 0 : whichStokes(0)=Stokes::I;
1860 0 : if (polType(0)=="X" || polType(0)=="Y") {
1861 0 : whichStokes(1)=Stokes::Q;
1862 : }
1863 : else {
1864 0 : whichStokes(1)=Stokes::V;
1865 : }
1866 0 : break;
1867 0 : case 3:
1868 0 : whichStokes.resize(3);
1869 0 : whichStokes(0)=Stokes::I;
1870 0 : whichStokes(1)=Stokes::Q;
1871 0 : whichStokes(2)=Stokes::U;
1872 0 : break;
1873 0 : case 4:
1874 0 : whichStokes.resize(4);
1875 0 : whichStokes(0)=Stokes::I;
1876 0 : whichStokes(1)=Stokes::Q;
1877 0 : whichStokes(2)=Stokes::U;
1878 0 : whichStokes(3)=Stokes::V;
1879 0 : break;
1880 : };
1881 : }
1882 0 : }
1883 0 : StokesCoordinate myStokes(whichStokes);
1884 :
1885 : // Now set up coordinates for image. If the shape is
1886 : // 5 dimensional then we add another StokesCoordinate
1887 0 : CoordinateSystem coordInfo;
1888 0 : coordInfo.addCoordinate(myRaDec);
1889 0 : coordInfo.addCoordinate(myStokes);
1890 0 : if(shape.nelements()==5) coordInfo.addCoordinate(myStokes);
1891 0 : coordInfo.addCoordinate(mySpectral);
1892 0 : return coordInfo;
1893 :
1894 0 : }
1895 :
1896 : /* This function sets up the Stokes labelling for the polarization planes of
1897 : 'cImage', the complex image that stores gridded correlations.
1898 :
1899 : For example, if 'V' is asked for, we need to grid RR and LL (npol = 2 pol planes in cImage).
1900 :
1901 : Called from the first call of ImageStokesImageUtil::cImage().
1902 :
1903 : Input : shape = shape of cImage.
1904 : polRep : Circular or Linear - present in the data.
1905 : Input/Output : coords = Input polarization labelling is that of the target image ('V')
1906 : Output polarization labelling is for the correlations ('RR','LL').
1907 : whichStokes = empty on input (first time).
1908 : Fill in the data correlation labels ('RR', 'LL')
1909 :
1910 : */
1911 2356 : CoordinateSystem StokesImageUtil::CStokesCoord(//const IPosition& shape,
1912 : const CoordinateSystem& coord,
1913 : Vector<Int>& whichStokes,
1914 : StokesImageUtil::PolRep polRep) {
1915 : /*
1916 : Int nx=shape(0);
1917 : Int ny=shape(1);
1918 : Int npol=shape(2);
1919 : Int nchan=shape(3);
1920 : AlwaysAssert(nx>0, AipsError);
1921 : AlwaysAssert(ny>0, AipsError);
1922 : AlwaysAssert(npol>0, AipsError);
1923 : AlwaysAssert(nchan>0, AipsError);
1924 : */
1925 :
1926 2356 : Int directionIndex=coord.findCoordinate(Coordinate::DIRECTION);
1927 2356 : DirectionCoordinate directionCoord=coord.directionCoordinate(directionIndex);
1928 :
1929 2356 : Int spectralIndex=coord.findCoordinate(Coordinate::SPECTRAL);
1930 2356 : SpectralCoordinate spectralCoord=coord.spectralCoordinate(spectralIndex);
1931 :
1932 2356 : Int stokesIndex=coord.findCoordinate(Coordinate::STOKES);
1933 : StokesCoordinate
1934 2356 : stokesCoord=coord.stokesCoordinate(stokesIndex);
1935 :
1936 : /* STOKESDBG */ //cout << "Util : CStokesCoord - input - stokesCoord : " << stokesCoord.stokes() << endl;
1937 :
1938 : // Polarization: If the specified whichStokes are ok, we use them
1939 : // if(Int(whichStokes.nelements())!=npol) {
1940 : // whichStokes.resize(npol);
1941 : // whichStokes=0;
1942 2356 : changeLabelsStokesToCorrStokes(stokesCoord, polRep, whichStokes);
1943 : //}
1944 2356 : AlwaysAssert(whichStokes.nelements(), AipsError);
1945 2356 : StokesCoordinate stokesCoordOut(whichStokes);
1946 :
1947 : /* STOKESDBG */ //cout << "Util : CStokesCoord - output - stokesCoord : " << stokesCoordOut.stokes() << endl;
1948 :
1949 : // Now set up coordinates
1950 2356 : CoordinateSystem coordInfo;
1951 2356 : coordInfo.addCoordinate(directionCoord);
1952 2356 : coordInfo.addCoordinate(stokesCoordOut);
1953 2356 : coordInfo.addCoordinate(spectralCoord);
1954 4712 : return coordInfo;
1955 :
1956 2356 : }
1957 :
1958 : /*
1959 : Logic : Read the desired image stokes from "stokesCoord" : 'V'
1960 : Read the desired "npol" from the length of stokesCoord : '2'
1961 : For npol = 1,2,3,4, if any image-pol choice requires an explicit mapping
1962 : to correlations, resize 'whichStokes' accordingly, and fill in.
1963 : For example, 'V' is mapped to 'RR,LL' for Circular, and 'XY,YX' for Linear.
1964 : Finally, check that all correlations in 'whichStokes' are present in the data.
1965 :
1966 : Subsequent calls to this function -- which already have the input stokesCoord in the
1967 : correlation basis, will not do anything.
1968 :
1969 : Rules :
1970 : npol = 1 : I, XX,YY,XY,YX,RR,LL.RL,LR --- choose only one correlation
1971 : Q, U, V -- map to 2 correlations as required for Circular vs Linear
1972 : Note : "I" is a special-case for which ftmachines and gridders know how to use only one plane.
1973 : npol = 2 : RRLL,RLLR,XXYY,XYYX; IV,QU for Circular; IQ,UV for Linear -- map to 2 correlations
1974 : IV, QU for Linear; IQ,UV for Circular -- map to 4 correlations
1975 : npol = 3,4 : Choose all 4 correlations.
1976 :
1977 : 'whichStokes' contains the output
1978 : */
1979 3557 : void StokesImageUtil::changeLabelsStokesToCorrStokes(StokesCoordinate &stokesCoord,
1980 : StokesImageUtil::PolRep polRep, //Int npol,
1981 : Vector<Int>&whichStokes)
1982 : {
1983 3557 : Int inputNPol = (stokesCoord.stokes()).nelements();
1984 3557 : Vector<Int> svec = stokesCoord.stokes();
1985 3557 : AlwaysAssert( (inputNPol==1 || inputNPol==2 || inputNPol==3 || inputNPol), AipsError );
1986 :
1987 3557 : switch( inputNPol )
1988 : {
1989 3357 : case 1:
1990 3357 : whichStokes.resize(1);
1991 : // by default, set to what the input is. This will take care of RR,LL,RL,LR,XX,YY,XY,YX and I.
1992 3357 : whichStokes[0]=svec[0];
1993 : // for Q,U,V, two correlations needs to be gridded
1994 3357 : if(polRep==StokesImageUtil::CIRCULAR && Stokes::type(svec[0])==Stokes::Q)
1995 22 : {whichStokes.resize(2); whichStokes[0]=Stokes::RL; whichStokes[1]=Stokes::LR; }
1996 3335 : else if(polRep==StokesImageUtil::CIRCULAR && Stokes::type(svec[0])==Stokes::U)
1997 52 : {whichStokes.resize(2); whichStokes[0]=Stokes::RL; whichStokes[1]=Stokes::LR; }
1998 3283 : else if(polRep==StokesImageUtil::CIRCULAR && Stokes::type(svec[0])==Stokes::V)
1999 6 : {whichStokes.resize(2); whichStokes[0]=Stokes::RR; whichStokes[1]=Stokes::LL; }
2000 3277 : else if(polRep==StokesImageUtil::LINEAR && Stokes::type(svec[0])==Stokes::Q)
2001 0 : {whichStokes.resize(2); whichStokes[0]=Stokes::XX; whichStokes[1]=Stokes::YY; }
2002 3277 : else if(polRep==StokesImageUtil::LINEAR && Stokes::type(svec[0])==Stokes::U)
2003 0 : {whichStokes.resize(2); whichStokes[0]=Stokes::XY; whichStokes[1]=Stokes::YX; }
2004 3277 : else if(polRep==StokesImageUtil::LINEAR && Stokes::type(svec[0])==Stokes::V)
2005 0 : {whichStokes.resize(2); whichStokes[0]=Stokes::XY; whichStokes[1]=Stokes::YX; }
2006 3357 : break;
2007 43 : case 2:
2008 43 : whichStokes.resize(2);
2009 : // by default, set to what the input is. This will take care of RRLL, RLLR, XXYY, XYYX
2010 43 : whichStokes(0)=svec[0];
2011 43 : whichStokes(1)=svec[1];
2012 : // only 2 correlations need to be gridded
2013 43 : if(polRep==StokesImageUtil::CIRCULAR && Stokes::type(svec[0])==Stokes::I && Stokes::type(svec[1])==Stokes::V)
2014 6 : {whichStokes(0)=Stokes::RR; whichStokes(1)=Stokes::LL;}
2015 37 : else if(polRep==StokesImageUtil::CIRCULAR && Stokes::type(svec[0])==Stokes::Q && Stokes::type(svec[1])==Stokes::U)
2016 6 : {whichStokes(0)=Stokes::RL; whichStokes(1)=Stokes::LR;}
2017 31 : else if(polRep==StokesImageUtil::LINEAR && Stokes::type(svec[0])==Stokes::I && Stokes::type(svec[1])==Stokes::Q)
2018 0 : {whichStokes(0)=Stokes::XX; whichStokes(1)=Stokes::YY;}
2019 31 : else if(polRep==StokesImageUtil::LINEAR && Stokes::type(svec[0])==Stokes::U && Stokes::type(svec[1])==Stokes::V)
2020 0 : {whichStokes(0)=Stokes::XY; whichStokes(1)=Stokes::YX;}
2021 : // all 4 correlations to be gridded. Only difference with above is CIRC to LIN
2022 37 : else if( (polRep==StokesImageUtil::LINEAR && Stokes::type(svec[0])==Stokes::I && Stokes::type(svec[1])==Stokes::V) ||
2023 6 : (polRep==StokesImageUtil::LINEAR && Stokes::type(svec[0])==Stokes::Q && Stokes::type(svec[1])==Stokes::U) )
2024 : {
2025 0 : whichStokes.resize(4);
2026 0 : whichStokes(0)=Stokes::XX; whichStokes(1)=Stokes::XY;
2027 0 : whichStokes(2)=Stokes::YX; whichStokes(3)=Stokes::YY;
2028 : }
2029 56 : else if( (polRep==StokesImageUtil::CIRCULAR && Stokes::type(svec[0])==Stokes::I && Stokes::type(svec[1])==Stokes::Q) ||
2030 25 : (polRep==StokesImageUtil::CIRCULAR && Stokes::type(svec[0])==Stokes::U && Stokes::type(svec[1])==Stokes::V) )
2031 : {
2032 0 : whichStokes.resize(4);
2033 0 : whichStokes(0)=Stokes::RR; whichStokes(1)=Stokes::RL;
2034 0 : whichStokes(2)=Stokes::LR; whichStokes(3)=Stokes::LL;
2035 : }
2036 43 : break;
2037 157 : case 3: /* If npol=3, all 4 correlations need to be gridded */
2038 : case 4: /* Select all 4 correlations when npol=4 */
2039 157 : whichStokes.resize(4);
2040 157 : if( polRep==StokesImageUtil::LINEAR)
2041 : {
2042 103 : whichStokes(0)=Stokes::XX;
2043 103 : whichStokes(1)=Stokes::XY;
2044 103 : whichStokes(2)=Stokes::YX;
2045 103 : whichStokes(3)=Stokes::YY;
2046 : }
2047 : else
2048 : {
2049 54 : whichStokes(0)=Stokes::RR;
2050 54 : whichStokes(1)=Stokes::RL;
2051 54 : whichStokes(2)=Stokes::LR;
2052 54 : whichStokes(3)=Stokes::LL;
2053 : }
2054 157 : break;
2055 : }
2056 :
2057 : // Verify that all entries in whichStokes are present in the data
2058 : // This is to catch things like : dataset contains only RR,LL, but the user has asked for IQUV.
2059 :
2060 :
2061 3557 : }// end of changeStokesToCorrStokes
2062 :
2063 :
2064 : /* This function is not called from anywhere !!!! It should go. */
2065 : #if 0
2066 : CoordinateSystem
2067 : StokesImageUtil::CStokesCoordFromImage(const ImageInterface<Complex>& image,
2068 : Vector<Int>& whichStokes,
2069 : StokesImageUtil::PolRep polRep) {
2070 :
2071 : IPosition shape=image.shape();
2072 :
2073 : Int npol=shape(2);
2074 :
2075 : CoordinateSystem coord=image.coordinates();
2076 :
2077 : Int directionIndex=coord.findCoordinate(Coordinate::DIRECTION);
2078 : DirectionCoordinate directionCoord=coord.directionCoordinate(directionIndex);
2079 :
2080 : Int spectralIndex=coord.findCoordinate(Coordinate::SPECTRAL);
2081 : SpectralCoordinate spectralCoord=coord.spectralCoordinate(spectralIndex);
2082 :
2083 : // Polarization: If the specified whichStokes are ok, we use them
2084 : // otherwise we guess
2085 : if(Int(whichStokes.nelements())!=npol) {
2086 : whichStokes.resize(npol);
2087 : // Polarization
2088 : if (polRep==StokesImageUtil::LINEAR) {
2089 : switch(npol) {
2090 : case 1:
2091 : whichStokes.resize(1);
2092 : whichStokes(1)=Stokes::I;
2093 : break;
2094 : case 2:
2095 : whichStokes.resize(2);
2096 : whichStokes(1)=Stokes::XX;
2097 : whichStokes(0)=Stokes::YY;
2098 : break;
2099 : default:
2100 : whichStokes.resize(4);
2101 : whichStokes(0)=Stokes::XX;
2102 : whichStokes(1)=Stokes::XY;
2103 : whichStokes(2)=Stokes::YX;
2104 : whichStokes(3)=Stokes::YY;
2105 : }
2106 : }
2107 : else {
2108 : switch(npol) {
2109 : case 1:
2110 : whichStokes.resize(1);
2111 : whichStokes(0)=Stokes::I;
2112 : break;
2113 : case 2:
2114 : whichStokes.resize(2);
2115 : whichStokes(0)=Stokes::LL;
2116 : whichStokes(1)=Stokes::RR;
2117 : break;
2118 : default:
2119 : whichStokes.resize(4);
2120 : whichStokes(0)=Stokes::LL;
2121 : whichStokes(1)=Stokes::LR;
2122 : whichStokes(2)=Stokes::RL;
2123 : whichStokes(3)=Stokes::RR;
2124 : };
2125 : }
2126 : }
2127 : StokesCoordinate stokesCoord(whichStokes);
2128 :
2129 : // Now set up coordinates for image.
2130 : CoordinateSystem coordInfo;
2131 : coordInfo.addCoordinate(directionCoord);
2132 : coordInfo.addCoordinate(stokesCoord);
2133 : coordInfo.addCoordinate(spectralCoord);
2134 : return coordInfo;
2135 :
2136 : }
2137 : #endif
2138 :
2139 : // Return a map to RA, DEC, pol, chan. 0 is RA, 1 is Dec, 2 is
2140 : // polarization, and 3 is Channel
2141 24338 : Bool StokesImageUtil::StokesMap(Vector<Int>& map, const CoordinateSystem& coord) {
2142 :
2143 24338 : map.resize(4);
2144 24338 : map=-1;
2145 :
2146 24338 : Int dirIndex=coord.findCoordinate(Coordinate::DIRECTION);
2147 24338 : if(dirIndex>-1) {
2148 24338 : Vector<Int> dirAxes=coord.pixelAxes(dirIndex);
2149 24338 : if(dirAxes.nelements()>2) {
2150 0 : return false;
2151 : }
2152 24338 : map(0)=dirAxes(0);
2153 24338 : map(1)=dirAxes(1);
2154 24338 : }
2155 : else {
2156 0 : return false;
2157 : }
2158 :
2159 24338 : Int stokesIndex=coord.findCoordinate(Coordinate::STOKES);
2160 24338 : if(stokesIndex>-1) {
2161 24338 : Vector<Int> stokesAxes=coord.pixelAxes(stokesIndex);
2162 24338 : if(stokesAxes.nelements()>1) {
2163 0 : return false;
2164 : }
2165 24338 : map(2)=stokesAxes(0);
2166 24338 : }
2167 : else {
2168 0 : return false;
2169 : }
2170 :
2171 24338 : Int spectralIndex=coord.findCoordinate(Coordinate::SPECTRAL);
2172 24338 : if(spectralIndex>-1) {
2173 24338 : Vector<Int> spectralAxes=coord.pixelAxes(spectralIndex);
2174 24338 : if(spectralAxes.nelements()>1) {
2175 0 : return false;
2176 : }
2177 24338 : map(3)=spectralAxes(0);
2178 24338 : }
2179 : else {
2180 0 : return false;
2181 : }
2182 :
2183 24338 : return true;
2184 : }
2185 :
2186 0 : void StokesImageUtil::BoxMask(ImageInterface<Float>& mask,
2187 : const IPosition& blc,
2188 : const IPosition& trc, const Float value)
2189 : {
2190 0 : LatticeStepper ls(mask.shape(), IPosition(4, mask.shape()(0), 1, 1, 1),
2191 0 : IPosition(4, 0, 1, 2, 3));
2192 0 : ls.subSection(blc, trc);
2193 0 : LatticeIterator<Float> mli(mask, ls);
2194 :
2195 : // Loop over all planes
2196 0 : for (mli.reset();!mli.atEnd();mli++) {
2197 0 : mli.rwCursor()=value;
2198 : }
2199 :
2200 0 : }
2201 :
2202 : // Change the representation used. The contents of the image
2203 : // are not changed!
2204 1201 : void StokesImageUtil::changeCStokesRep(ImageInterface<Complex>& image,
2205 : StokesImageUtil::PolRep polRep) {
2206 :
2207 1201 : CoordinateSystem coords=image.coordinates();
2208 :
2209 1201 : Int stokesIndex=coords.findCoordinate(Coordinate::STOKES);
2210 : StokesCoordinate
2211 1201 : stokesCoord=coords.stokesCoordinate(stokesIndex);
2212 :
2213 : ///// Int npol=stokesCoord.stokes().nelements(); // apparently not used and no side effects
2214 :
2215 : /* STOKESDBG */ //cout << "Util::changeCStokesRep - input - stokescoord : " << stokesCoord.stokes() << " npol : " << npol << endl;
2216 :
2217 1201 : Vector<Int> whichStokes(0);
2218 : //whichStokes=0;
2219 1201 : changeLabelsStokesToCorrStokes(stokesCoord, polRep, whichStokes);
2220 :
2221 1201 : StokesCoordinate newStokesCoord(whichStokes);
2222 1201 : coords.replaceCoordinate(newStokesCoord, stokesIndex);
2223 1201 : image.setCoordinateInfo(coords);
2224 : /* STOKESDBG */ //cout << "Util::changeCStokesRep - output - stokescoord : " << newStokesCoord.stokes() << endl;
2225 1201 : }
2226 :
2227 0 : Bool StokesImageUtil::
2228 : standardImageCoordinates(const ImageInterface<Complex>& image) {
2229 0 : return (standardImageCoordinates(image.coordinates()));
2230 : };
2231 :
2232 0 : Bool StokesImageUtil::
2233 : standardImageCoordinates(const ImageInterface<Float>& image) {
2234 0 : return (standardImageCoordinates(image.coordinates()));
2235 : };
2236 :
2237 0 : Bool StokesImageUtil::
2238 : standardImageCoordinates(const CoordinateSystem& coords) {
2239 0 : Bool isStandard = true;
2240 : {
2241 0 : Int ind=coords.findCoordinate(Coordinate::DIRECTION);
2242 0 : if (ind != 0) isStandard = false;
2243 0 : ind=coords.findCoordinate(Coordinate::STOKES);
2244 0 : if (ind != 1) isStandard = false;
2245 0 : ind=coords.findCoordinate(Coordinate::SPECTRAL);
2246 0 : if (ind != 2) isStandard = false;
2247 : }
2248 0 : return isStandard;
2249 : };
2250 :
2251 : } //#End casa namespace
2252 :
2253 :
|
