Line data Source code
1 : // -*- C++ -*-
2 : //# Utils.cc: Implementation of global functions from Utils.h
3 : //# Copyright (C) 1997,1998,1999,2000,2001,2002,2003
4 : //# Associated Universities, Inc. Washington DC, USA.
5 : //#
6 : //# This library is free software; you can redistribute it and/or modify it
7 : //# under the terms of the GNU Library General Public License as published by
8 : //# the Free Software Foundation; either version 2 of the License, or (at your
9 : //# option) any later version.
10 : //#
11 : //# This library is distributed in the hope that it will be useful, but WITHOUT
12 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
13 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
14 : //# License for more details.
15 : //#
16 : //# You should have received a copy of the GNU Library General Public License
17 : //# along with this library; if not, write to the Free Software Foundation,
18 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
19 : //#
20 : //# Correspondence concerning AIPS++ should be addressed as follows:
21 : //# Internet email: casa-feedback@nrao.edu.
22 : //# Postal address: AIPS++ Project Office
23 : //# National Radio Astronomy Observatory
24 : //# 520 Edgemont Road
25 : //# Charlottesville, VA 22903-2475 USA
26 : //#
27 : //# $Id$
28 : #include <casacore/ms/MeasurementSets/MSRange.h>
29 : #include <msvis/MSVis/VisBuffer2.h>
30 : #include <casacore/casa/Logging/LogIO.h>
31 : #include <casacore/ms/MeasurementSets/MSColumns.h>
32 : #include <casacore/measures/Measures/MEpoch.h>
33 : #include <casacore/measures/Measures/MeasTable.h>
34 : #include <synthesis/TransformMachines2/Utils.h>
35 : #include <synthesis/TransformMachines/StokesImageUtil.h>
36 : #include <synthesis/Utilities/FFT2D.h>
37 : #include <casacore/casa/Utilities/Assert.h>
38 : #include <casacore/casa/Arrays/Vector.h>
39 : #include <casacore/casa/Arrays/ArrayMath.h>
40 : #include <casacore/lattices/LEL/LatticeExpr.h>
41 : #include <casacore/images/Images/PagedImage.h>
42 : #include <casacore/images/Images/ImageRegrid.h>
43 : #include <casacore/casa/Containers/Record.h>
44 : #include <casacore/lattices/Lattices/LatticeIterator.h>
45 : #include <casacore/lattices/Lattices/TiledLineStepper.h>
46 : #include <casacore/lattices/Lattices/LatticeStepper.h>
47 : #include <casacore/lattices/LatticeMath/LatticeFFT.h>
48 : #include <casacore/casa/OS/Timer.h>
49 :
50 : #include <casacore/casa/System/Aipsrc.h>
51 : #include <msvis/MSVis/VisibilityIterator2.h>
52 :
53 : using namespace casacore;
54 : namespace casa{
55 : using namespace vi;
56 : namespace refim{
57 : //
58 : //--------------------------------------------------------------------------------------------
59 : //
60 0 : void storeImg(String fileName,ImageInterface<Complex>& theImg, Bool writeReIm)
61 : {
62 0 : PagedImage<Complex> ctmp(theImg.shape(), theImg.coordinates(), fileName);
63 0 : LatticeExpr<Complex> le(theImg);
64 0 : ctmp.copyData(le);
65 0 : if (writeReIm)
66 : {
67 0 : ostringstream reName,imName;
68 0 : reName << "re" << fileName;
69 0 : imName << "im" << fileName;
70 : {
71 0 : PagedImage<Float> tmp(theImg.shape(), theImg.coordinates(), reName);
72 : //LatticeExpr<Float> le(abs(theImg));
73 0 : LatticeExpr<Float> le(real(theImg));
74 0 : tmp.copyData(le);
75 0 : }
76 : {
77 0 : PagedImage<Float> tmp(theImg.shape(), theImg.coordinates(), imName);
78 0 : LatticeExpr<Float> le(arg(theImg));
79 0 : tmp.copyData(le);
80 0 : }
81 0 : }
82 0 : }
83 :
84 0 : void storeImg(String fileName,ImageInterface<Float>& theImg)
85 : {
86 0 : PagedImage<Float> tmp(theImg.shape(), theImg.coordinates(), fileName);
87 0 : LatticeExpr<Float> le(theImg);
88 0 : tmp.copyData(le);
89 0 : }
90 :
91 0 : void storeArrayAsImage(String fileName, const CoordinateSystem& coord,
92 : const Array<Complex>& theImg)
93 : {
94 0 : PagedImage<Complex> ctmp(theImg.shape(), coord, fileName);
95 0 : ctmp.put(theImg);
96 0 : }
97 0 : void storeArrayAsImage(String fileName, const CoordinateSystem& coord,
98 : const Array<DComplex>& theImg)
99 : {
100 0 : PagedImage<DComplex> ctmp(theImg.shape(), coord, fileName);
101 0 : ctmp.put(theImg);
102 0 : }
103 0 : void storeArrayAsImage(String fileName, const CoordinateSystem& coord,
104 : const Array<Float>& theImg)
105 : {
106 0 : PagedImage<Float> ctmp(theImg.shape(), coord, fileName);
107 0 : ctmp.put(theImg);
108 0 : }
109 : //
110 : //---------------------------------------------------------------------
111 : //
112 : // Get the time stamp for the for the current visibility integration.
113 : // Since VisTimeAverager() does not accumulate auto-correlations (it
114 : // should though!), and the VisBuffer can potentially have
115 : // auto-correlation placeholders, vb.time()(0) may not be correct (it
116 : // will be in fact zero when AC rows are present). So look for the
117 : // first timestamp of a row corresponding to an unflagged
118 : // cross-correlation.
119 : //
120 151332 : Double getCurrentTimeStamp(const VisBuffer2& vb)
121 : {
122 : //LogIO os(LogOrigin("Utils", "getCurrentTimeStamp", WHERE));
123 :
124 151332 : Int N=vb.nRows();
125 172584 : for(Int i=0;i<N;i++)
126 : {
127 172584 : if ((!vb.flagRow()(i)) && (vb.antenna1()(i) != vb.antenna2()(i)))
128 151332 : return vb.time()(i);
129 : }
130 : //os << "No unflagged row found! This is a major problem/bug" << LogIO::WARN;
131 0 : return 0.0;
132 : }
133 : //
134 : //---------------------------------------------------------------------
135 : // Compute the Parallactic Angle for the give VisBuffer
136 : //
137 0 : void getHADec(MeasurementSet& ms, const VisBuffer2& vb,
138 : Double& HA, Double& RA, Double& Dec)
139 : {
140 0 : MEpoch last;
141 0 : Double time = getCurrentTimeStamp(vb);
142 :
143 0 : Unit Second("s"), Day("d");
144 0 : MEpoch epoch(Quantity(time,Second),MEpoch::TAI);
145 0 : MPosition pos;
146 0 : MSObservationColumns msoc(ms.observation());
147 0 : String ObsName=msoc.telescopeName()(vb.arrayId()(0));
148 :
149 0 : if (!MeasTable::Observatory(pos,ObsName))
150 0 : throw(AipsError("Observatory position for "+ ObsName + " not found"));
151 :
152 0 : MeasFrame frame(epoch,pos);
153 0 : MEpoch::Convert toLAST = MEpoch::Convert(MEpoch::Ref(MEpoch::TAI,frame),
154 0 : MEpoch::Ref(MEpoch::LAST,frame));
155 0 : RA=vb.direction1()(0).getAngle().getValue()(0);
156 0 : if (RA < 0.0) RA += M_PI*2.0;
157 0 : Dec=vb.direction1()(0).getAngle().getValue()(1);
158 :
159 0 : last = toLAST(epoch);
160 0 : Double LST = last.get(Day).getValue();
161 0 : LST -= floor(LST); // Extract the fractional day
162 0 : LST *= 2*C::pi;// Convert to Raidan
163 :
164 0 : cout << "LST = " << LST << " " << LST/(2*C::pi) << endl;
165 :
166 0 : HA = LST - RA;
167 0 : }
168 : //
169 : //---------------------------------------------------------------------
170 : // Compute the Parallactic Angle for the give VisBuffer
171 : //
172 151332 : Double getPA(const vi::VisBuffer2& vb)
173 : {
174 151332 : return (Double)(vb.feedPa(getCurrentTimeStamp(vb))(0));
175 : // Double pa=0;
176 : // Int n=0;
177 : // Vector<Float> antPA = vb.feed_pa(getCurrentTimeStamp(vb));
178 : // for (uInt i=0;i<antPA.nelements();i++)
179 : // {
180 : // if (!vb.msColumns().antenna().flagRow()(i))
181 : // {pa += antPA(i); n++;break;}
182 : // }
183 : // // pa = sum(antPA)/(antPA.nelements()-1);
184 : // if (n==0)
185 : // throw(AipsError("No unflagged antenna found in getPA()."));
186 : // return pa/n;
187 : }
188 : //
189 : //---------------------------------------------------------------------
190 : //
191 : // Make stokes axis, given the polarization types.
192 : //
193 0 : void makeStokesAxis(Int npol_p, Vector<String>& polType, Vector<Int>& whichStokes)
194 : {
195 : // Vector<String> polType=msc.feed().polarizationType()(0);
196 : StokesImageUtil::PolRep polRep_p;
197 0 : LogIO os(LogOrigin("Utils", "makeStokesAxis", WHERE));
198 :
199 0 : if (polType(0)!="X" && polType(0)!="Y" &&
200 0 : polType(0)!="R" && polType(0)!="L")
201 : {
202 : os << "Warning: Unknown stokes types in feed table: ["
203 0 : << polType(0) << ", " << polType(1) << "]" << endl
204 0 : << "Results open to question!" << LogIO::POST;
205 : }
206 :
207 0 : if (polType(0)=="X" || polType(0)=="Y")
208 : {
209 0 : polRep_p=StokesImageUtil::LINEAR;
210 0 : os << "Preferred polarization representation is linear" << LogIO::POST;
211 : }
212 : else
213 : {
214 0 : polRep_p=StokesImageUtil::CIRCULAR;
215 0 : os << "Preferred polarization representation is circular" << LogIO::POST;
216 : }
217 :
218 : // Vector<Int> whichStokes(npol_p);
219 0 : switch(npol_p)
220 : {
221 0 : case 1:
222 0 : whichStokes.resize(1);
223 0 : whichStokes(0)=Stokes::I;
224 0 : os << "Image polarization = Stokes I" << LogIO::POST;
225 0 : break;
226 0 : case 2:
227 0 : whichStokes.resize(2);
228 0 : whichStokes(0)=Stokes::I;
229 0 : if (polRep_p==StokesImageUtil::LINEAR)
230 : {
231 0 : whichStokes(1)=Stokes::Q;
232 0 : os << "Image polarization = Stokes I,Q" << LogIO::POST;
233 : }
234 : else
235 : {
236 0 : whichStokes(1)=Stokes::V;
237 0 : os << "Image polarization = Stokes I,V" << LogIO::POST;
238 : }
239 0 : break;
240 0 : case 3:
241 0 : whichStokes.resize(3);
242 0 : whichStokes(0)=Stokes::I;
243 0 : whichStokes(1)=Stokes::Q;
244 0 : whichStokes(1)=Stokes::U;
245 0 : os << "Image polarization = Stokes I,Q,U" << LogIO::POST;
246 0 : break;
247 0 : case 4:
248 0 : whichStokes.resize(4);
249 0 : whichStokes(0)=Stokes::I;
250 0 : whichStokes(1)=Stokes::Q;
251 0 : whichStokes(2)=Stokes::U;
252 0 : whichStokes(3)=Stokes::V;
253 0 : os << "Image polarization = Stokes I,Q,U,V" << LogIO::POST;
254 0 : break;
255 0 : default:
256 : os << LogIO::SEVERE << "Illegal number of Stokes parameters: " << npol_p
257 0 : << LogIO::POST;
258 : };
259 0 : }
260 : //
261 : //--------------------------------------------------------------------------------------------
262 : //
263 0 : Bool isVBNaN(const VisBuffer2 &vb,String& mesg)
264 : {
265 0 : IPosition ndx(3,0);
266 0 : Int N = vb.nRows();
267 0 : for(ndx(2)=0;ndx(2)<N;ndx(2)++)
268 0 : if (std::isnan(vb.visCubeModel()(ndx).real()) ||
269 0 : std::isnan(vb.visCubeModel()(ndx).imag())
270 : )
271 : {
272 0 : ostringstream os;
273 0 : os << ndx(2) << " " << vb.antenna1()(ndx(2)) << "-" << vb.antenna2()(ndx(2)) << " "
274 0 : << vb.flagCube()(ndx) << " "
275 : //<< vb.flag()(0,ndx(2)) << " "
276 0 : << vb.weight();
277 0 : mesg = os.str().c_str();
278 0 : return true;
279 0 : }
280 0 : return false;
281 0 : }
282 : //
283 : //--------------------------------------------------------------------------------------------
284 : //
285 : /////////////////////////////////////////////////////////////////////////////
286 : //
287 : // IChangeDetector - an interface class to detect changes in the VisBuffer
288 : // Exact meaning of the "change" is defined in the derived classes
289 : // (e.g. a change in the parallactic angle)
290 :
291 : // return true if a change occurs somewhere in the buffer
292 : using namespace vi;
293 0 : Bool IChangeDetector::changed(const VisBuffer2 &vb) const
294 : {
295 0 : for (rownr_t i=0;i<vb.nRows();++i)
296 0 : if (changed(vb,i)) return true;
297 0 : return false;
298 : }
299 :
300 : // return true if a change occurs somewhere in the buffer starting from row1
301 : // up to row2 (row2=-1 means up to the end of the buffer). The row number,
302 : // where the change occurs is returned in the row2 parameter
303 0 : Bool IChangeDetector::changedBuffer(const VisBuffer2 &vb, Int row1,
304 : Int &row2) const
305 : {
306 0 : if (row1<0) row1=0;
307 0 : int jrow = row2;
308 0 : if (jrow < 0) jrow = vb.nRows()-1;
309 0 : DebugAssert(jrow<(int)vb.nRows(),AipsError);
310 :
311 : // It is not important now to have a separate function for a "block"
312 : // operation. Because an appropriate caching is implemented inside
313 : // VisBuffer, changed(vb,row) can be called in the
314 : // loop without a perfomance penalty. This method returns the
315 : // first row where the change occured rather than the last unchanged
316 : // row as it was for BeamSkyJones::changedBuffer
317 :
318 0 : for (int ii=row1;ii<=jrow;++ii)
319 0 : if (changed(vb,ii)) {
320 0 : row2 = ii;
321 0 : return true;
322 : }
323 0 : return false;
324 : }
325 :
326 : // a virtual destructor to make the compiler happy
327 352 : IChangeDetector::~IChangeDetector() {}
328 :
329 : //
330 : /////////////////////////////////////////////////////////////////////////////
331 :
332 : /////////////////////////////////////////////////////////////////////////////
333 : //
334 : // ParAngleChangeDetector - a class to detect a change in the parallatic
335 : // angle
336 : //
337 :
338 : // set up the tolerance, which determines how much the position angle should
339 : // change to report the change by this class
340 0 : ParAngleChangeDetector::ParAngleChangeDetector(const Quantity &pa_tolerance)
341 0 : : pa_tolerance_p(pa_tolerance.getValue("rad")),
342 0 : last_pa_p(1000.) {} // 1000 is >> 2pi, so it is changed
343 : // after construction
344 :
345 : // Set the value of the PA tolerance
346 93 : void ParAngleChangeDetector::setTolerance(const Quantity &pa_tolerance)
347 : {
348 93 : pa_tolerance_p = (pa_tolerance.getValue("rad"));
349 93 : }
350 : // reset to the state which exist just after construction
351 566 : void ParAngleChangeDetector::reset()
352 : {
353 566 : last_pa_p=1000.; // it is >> 2pi, which would force a changed state
354 566 : }
355 :
356 : // return parallactic angle tolerance
357 34607 : Quantity ParAngleChangeDetector::getParAngleTolerance() const
358 : {
359 34607 : return Quantity(pa_tolerance_p,"rad");
360 : }
361 :
362 : // return true if a change occurs in the given row since the last call
363 : // of update
364 27989 : Bool ParAngleChangeDetector::changed(const VisBuffer2 &vb, Int row) const
365 : {
366 27989 : if (row<0) row=0;
367 : // const Double feed1_pa=vb.feed1_pa()[row];
368 27989 : Double feed1_pa=getPA(vb);
369 27989 : if (abs(feed1_pa-last_pa_p) > pa_tolerance_p)
370 : {
371 : // cout << "Utils: " << feed1_pa*57.295 << " " << last_pa_p*57.295 << " " << abs(feed1_pa-last_pa_p)*57.295 << " " << ttt*57.295 << " " << vb.time()(0)-4.51738e+09 << endl;
372 422 : return true;
373 : }
374 : // const Double feed2_pa=vb.feed2_pa()[row];
375 27567 : Double feed2_pa = getPA(vb);
376 27567 : if (abs(feed2_pa-last_pa_p) > pa_tolerance_p)
377 : {
378 : // cout << "Utils: " << feed2_pa*57.295 << " "
379 : // << last_pa_p*57.295 << " "
380 : // << abs(feed2_pa-last_pa_p)*57.295 << " " << ttt*57.295 << vb.time()(0)-4.51738e+09 <<endl;
381 0 : return true;
382 : }
383 27567 : return false;
384 : }
385 :
386 : // start looking for a change from the given row of the VisBuffer
387 211 : void ParAngleChangeDetector::update(const VisBuffer2 &vb, Int row)
388 : {
389 211 : if (row<0) row=0;
390 211 : const Double feed1_pa=vb.feedPa1()[row];
391 211 : const Double feed2_pa=vb.feedPa2()[row];
392 211 : if (abs(feed1_pa-feed2_pa)>pa_tolerance_p) {
393 0 : LogIO os;
394 0 : os<<LogIO::WARN << LogOrigin("ParAngleChangeDetector","update")
395 : <<"The parallactic angle is different at different stations"
396 0 : <<LogIO::POST<<LogIO::WARN <<LogOrigin("ParAngleChangeDetector","update")
397 0 : <<"The result may be incorrect. Continuing anyway."<<LogIO::POST;
398 0 : }
399 211 : last_pa_p=(feed1_pa+feed2_pa)/2.;
400 211 : }
401 :
402 0 : Int getPhaseCenter(MeasurementSet& ms, MDirection& dir0, Int whichField)
403 : {
404 0 : MSFieldColumns msfc(ms.field());
405 0 : if (whichField < 0)
406 : {
407 0 : MSRange msr(ms);
408 : //
409 : // An array of shape [2,1,1]!
410 : //
411 0 : Vector<Int> fieldId;
412 0 : fieldId = msr.range(MSS::FIELD_ID).asArrayInt(RecordFieldId(0));
413 :
414 0 : Array<Double> phaseDir = msfc.phaseDir().getColumn();
415 :
416 0 : if (fieldId.nelements() <= 0)
417 0 : throw(AipsError("getPhaseCenter: No fields found in the selected MS."));
418 :
419 0 : IPosition ndx0(3,0,0,0),ndx1(3,1,0,0);
420 :
421 : Double maxRA, maxDec, RA,Dec,RA0,Dec0, dist0;
422 0 : maxRA = maxDec=0;
423 0 : for(uInt i=0;i<fieldId.nelements();i++)
424 : {
425 0 : RA = phaseDir(IPosition(3,0,0,fieldId(i)));
426 0 : Dec = phaseDir(IPosition(3,1,0,fieldId(i)));
427 0 : maxRA += RA; maxDec += Dec;
428 : }
429 0 : RA0=maxRA/fieldId.nelements(); Dec0=maxDec/fieldId.nelements();
430 :
431 0 : dist0=10000;
432 0 : Int field0=0;
433 0 : for(uInt i=0;i<fieldId.nelements();i++)
434 : {
435 0 : RA = RA0-phaseDir(IPosition(3,0,0,fieldId(i)));
436 0 : Dec = Dec0-phaseDir(IPosition(3,1,0,fieldId(i)));
437 0 : Double dist=sqrt(RA*RA + Dec*Dec);
438 0 : if (dist < dist0) {field0=fieldId(i);dist0=dist;};
439 : }
440 0 : dir0=msfc.phaseDirMeasCol()(field0)(IPosition(1,0));
441 : //dir0=msfc.phaseDirMeasCol()(6)(IPosition(1,0));
442 0 : return field0;
443 0 : }
444 : else
445 : {
446 0 : dir0=msfc.phaseDirMeasCol()(whichField)(IPosition(1,0));
447 0 : return whichField;
448 : }
449 0 : }
450 : //
451 : //------------------------------------------------------------------
452 : //
453 0 : Bool findMaxAbsLattice(const ImageInterface<Float>& lattice,
454 : Float& maxAbs,IPosition& posMaxAbs)
455 : {
456 0 : posMaxAbs = IPosition(lattice.shape().nelements(), 0);
457 0 : maxAbs=0.0;
458 :
459 0 : const IPosition tileShape = lattice.niceCursorShape();
460 0 : TiledLineStepper ls(lattice.shape(), tileShape, 0);
461 : {
462 0 : RO_LatticeIterator<Float> li(lattice, ls);
463 0 : for(li.reset();!li.atEnd();li++)
464 : {
465 0 : IPosition posMax=li.position();
466 0 : IPosition posMin=li.position();
467 0 : Float maxVal=0.0;
468 0 : Float minVal=0.0;
469 :
470 0 : minMax(minVal, maxVal, posMin, posMax, li.cursor());
471 0 : if((maxVal)>(maxAbs))
472 : {
473 0 : maxAbs=maxVal;
474 0 : posMaxAbs=li.position();
475 0 : posMaxAbs(0)=posMax(0);
476 : }
477 0 : }
478 0 : }
479 :
480 0 : return true;
481 0 : };
482 : //
483 : //------------------------------------------------------------------
484 : //
485 0 : Bool findMaxAbsLattice(const ImageInterface<Float>& masklat,
486 : const Lattice<Float>& lattice,
487 : Float& maxAbs,IPosition& posMaxAbs,
488 : Bool flip)
489 : {
490 :
491 0 : AlwaysAssert(masklat.shape()==lattice.shape(), AipsError);
492 :
493 0 : Array<Float> msk;
494 :
495 0 : posMaxAbs = IPosition(lattice.shape().nelements(), 0);
496 0 : maxAbs=0.0;
497 : //maxAbs=-1.0e+10;
498 0 : const IPosition tileShape = lattice.niceCursorShape();
499 0 : TiledLineStepper ls(lattice.shape(), tileShape, 0);
500 0 : TiledLineStepper lsm(masklat.shape(), tileShape, 0);
501 : {
502 0 : RO_LatticeIterator<Float> li(lattice, ls);
503 0 : RO_LatticeIterator<Float> lim(masklat, lsm);
504 0 : for(li.reset(),lim.reset();!li.atEnd();li++,lim++)
505 : {
506 0 : IPosition posMax=li.position();
507 0 : IPosition posMin=li.position();
508 0 : Float maxVal=0.0;
509 0 : Float minVal=0.0;
510 :
511 0 : msk = lim.cursor();
512 0 : if(flip) msk = (Float)1.0 - msk;
513 :
514 : //minMaxMasked(minVal, maxVal, posMin, posMax, li.cursor(),lim.cursor());
515 0 : minMaxMasked(minVal, maxVal, posMin, posMax, li.cursor(),msk);
516 :
517 :
518 0 : if((maxVal)>(maxAbs))
519 : {
520 0 : maxAbs=maxVal;
521 0 : posMaxAbs=li.position();
522 0 : posMaxAbs(0)=posMax(0);
523 : }
524 0 : }
525 0 : }
526 :
527 0 : return true;
528 0 : }
529 : //
530 :
531 :
532 : //---------------------------------------------------------------
533 : //Rotate a complex array using a the given coordinate system and the
534 : //angle in radians. Default interpolation method is "CUBIC".
535 : //Axeses corresponding to Linear coordinates in the given
536 : //CoordinateSystem object are rotated. Rotation is done using
537 : //ImageRegrid object, about the pixel given by (N+1)/2 where N is
538 : //the number of pixels along the axis.
539 : //
540 0 : void SynthesisUtils::rotateComplexArray(LogIO& logio, Array<Complex>& inArray,
541 : CoordinateSystem& inCS,
542 : Array<Complex>& outArray,
543 : Double dAngleRad,
544 : String interpMethod,
545 : Bool modifyInCS)
546 : {
547 : // logio << LogOrigin("SynthesisUtils", "rotateComplexArray")
548 : // << "Rotating CF using " << interpMethod << " interpolation."
549 : // << LogIO::POST;
550 : (void)logio;
551 : //
552 : // If no rotation required, just copy the inArray to outArray.
553 : //
554 : // if (abs(dAngleRad) < 0.1)
555 :
556 : // IPosition tt;
557 : // inCS.list(logio,MDoppler::RADIO,tt,tt);
558 :
559 0 : if (abs(dAngleRad) == 0.0)
560 : {
561 0 : outArray.reference(inArray);
562 : // outArray.assign(inArray);
563 0 : return;
564 : }
565 : //
566 : // Re-grid inImage onto outImage
567 : //
568 0 : Vector<Int> pixelAxes;
569 0 : Int coordInd = -1;
570 : // Extract LINRAR coords from inCS.
571 : // Extract axes2
572 :
573 0 : if(modifyInCS){
574 0 : Vector<Double> refPix = inCS.referencePixel();
575 0 : refPix(0) = (Int)((inArray.shape()(0))/2.0);
576 0 : refPix(1) = (Int)((inArray.shape()(1))/2.0);
577 0 : inCS.setReferencePixel(refPix);
578 0 : }
579 :
580 0 : coordInd = inCS.findCoordinate(Coordinate::LINEAR);
581 0 : Bool haveLinear = true;
582 :
583 0 : if(coordInd == -1){ // no linear coordinate found, look for DIRECTION instead
584 0 : coordInd = inCS.findCoordinate(Coordinate::DIRECTION);
585 0 : haveLinear = false;
586 : }
587 :
588 0 : pixelAxes=inCS.pixelAxes(coordInd);
589 0 : IPosition axes2(pixelAxes);
590 : // Set linear transformation matrix in inCS.
591 : // CoordinateSystem outCS =
592 : // ImageRegrid<Complex>::makeCoordinateSystem (logio, outCS, inCS, axes2);
593 :
594 0 : CoordinateSystem outCS(inCS);
595 :
596 0 : Matrix<Double> xf = outCS.coordinate(coordInd).linearTransform();
597 0 : Matrix<Double> rotm(2,2);
598 0 : rotm(0,0) = cos(dAngleRad); rotm(0,1) = sin(dAngleRad);
599 0 : rotm(1,0) = -rotm(0,1); rotm(1,1) = rotm(0,0);
600 :
601 : // Double s = sin(dAngleRad);
602 : // Double c = cos(dAngleRad);
603 : // rotm(0,0) = c; rotm(0,1) = s;
604 : // rotm(1,0) = -s; rotm(1,1) = c;
605 :
606 : // Create new linear transform matrix
607 0 : Matrix<Double> xform(2,2);
608 0 : xform(0,0) = rotm(0,0)*xf(0,0)+rotm(0,1)*xf(1,0);
609 0 : xform(0,1) = rotm(0,0)*xf(0,1)+rotm(0,1)*xf(1,1);
610 0 : xform(1,0) = rotm(1,0)*xf(0,0)+rotm(1,1)*xf(1,0);
611 0 : xform(1,1) = rotm(1,0)*xf(0,1)+rotm(1,1)*xf(1,1);
612 :
613 0 : if(haveLinear){
614 0 : LinearCoordinate linCoords = outCS.linearCoordinate(coordInd);
615 0 : linCoords.setLinearTransform(xform);
616 0 : outCS.replaceCoordinate(linCoords, coordInd);
617 0 : }
618 : else{
619 0 : DirectionCoordinate dirCoords = outCS.directionCoordinate(coordInd);
620 0 : dirCoords.setLinearTransform(xform);
621 0 : outCS.replaceCoordinate(dirCoords, coordInd);
622 0 : }
623 :
624 0 : outArray.resize(inArray.shape());
625 0 : outArray.set(0);
626 : //
627 : // Make an image out of inArray and inCS --> inImage
628 : //
629 : // TempImage<Complex> inImage(inArray.shape(), inCS);
630 : {
631 0 : TempImage<Float> inImage(inArray.shape(),inCS);
632 0 : TempImage<Float> outImage(outArray.shape(), outCS);
633 0 : ImageRegrid<Float> ir;
634 0 : Interpolate2D::Method interpolationMethod = Interpolate2D::stringToMethod(interpMethod);
635 : //------------------------------------------------------------------------
636 : // Rotated the real part
637 : //
638 0 : inImage.copyData(LatticeExpr<Float>(real(ArrayLattice<Complex>(inArray))));
639 0 : outImage.set(0.0);
640 :
641 0 : ir.regrid(outImage, interpolationMethod, axes2, inImage);
642 0 : setReal(outArray,outImage.get());
643 : //------------------------------------------------------------------------
644 : // Rotated the imaginary part
645 : //
646 0 : inImage.copyData(LatticeExpr<Float>(imag(ArrayLattice<Complex>(inArray))));
647 0 : outImage.set(0.0);
648 :
649 0 : ir.regrid(outImage, interpolationMethod, axes2, inImage);
650 0 : setImag(outArray,outImage.get());
651 0 : }
652 0 : }
653 : //
654 : //---------------------------------------------------------------
655 : //
656 0 : void SynthesisUtils::findLatticeMax(const ImageInterface<Complex>& lattice,
657 : Vector<Float>& maxAbs,
658 : Vector<IPosition>& posMaxAbs)
659 : {
660 0 : IPosition lshape(lattice.shape());
661 0 : IPosition ndx(lshape);
662 0 : Int nPol=lshape(2);
663 0 : posMaxAbs.resize(nPol);
664 0 : for(Int i=0;i<nPol;i++)
665 0 : posMaxAbs(i)=IPosition(lattice.shape().nelements(), 0);
666 0 : maxAbs.resize(nPol);
667 0 : ndx=0;
668 :
669 0 : for(Int s2=0;s2<lshape(2);s2++)
670 0 : for(Int s3=0;s3<lshape(3);s3++)
671 : {
672 0 : ndx(2) = s2; ndx(3)=s3;
673 : {
674 : //
675 : // Locate the pixel with the peak value. That's the
676 : // origin in pixel co-ordinates.
677 : //
678 0 : maxAbs(s2)=0;
679 0 : posMaxAbs(s2) = 0;
680 0 : for(ndx(1)=0;ndx(1)<lshape(1);ndx(1)++)
681 0 : for(ndx(0)=0;ndx(0)<lshape(0);ndx(0)++)
682 0 : if (abs(lattice(ndx)) > maxAbs(s2))
683 0 : {posMaxAbs(s2) = ndx;maxAbs(s2)=abs(lattice(ndx));}
684 : }
685 : }
686 0 : }
687 : //
688 : //---------------------------------------------------------------
689 : //
690 0 : void SynthesisUtils::findLatticeMax(const Array<Complex>& lattice,
691 : Vector<Float>& maxAbs,
692 : Vector<IPosition>& posMaxAbs)
693 : {
694 0 : IPosition lshape(lattice.shape());
695 0 : IPosition ndx(lshape);
696 0 : Int nPol=lshape(2);
697 0 : posMaxAbs.resize(nPol);
698 0 : for(Int i=0;i<nPol;i++)
699 0 : posMaxAbs(i)=IPosition(lattice.shape().nelements(), 0);
700 0 : maxAbs.resize(nPol);
701 0 : ndx=0;
702 :
703 0 : for(Int s2=0;s2<lshape(2);s2++)
704 0 : for(Int s3=0;s3<lshape(3);s3++)
705 : {
706 0 : ndx(2) = s2; ndx(3)=s3;
707 : {
708 : //
709 : // Locate the pixel with the peak value. That's the
710 : // origin in pixel co-ordinates.
711 : //
712 0 : maxAbs(s2)=0;
713 0 : posMaxAbs(s2) = 0;
714 0 : for(ndx(1)=0;ndx(1)<lshape(1);ndx(1)++)
715 0 : for(ndx(0)=0;ndx(0)<lshape(0);ndx(0)++)
716 0 : if (abs(lattice(ndx)) > maxAbs(s2))
717 0 : {posMaxAbs(s2) = ndx;maxAbs(s2)=abs(lattice(ndx));}
718 : }
719 : }
720 0 : }
721 : //
722 : //---------------------------------------------------------------
723 : //
724 0 : void SynthesisUtils::findLatticeMax(const ImageInterface<Float>& lattice,
725 : Vector<Float>& maxAbs,
726 : Vector<IPosition>& posMaxAbs)
727 : {
728 0 : IPosition lshape(lattice.shape());
729 0 : IPosition ndx(lshape);
730 0 : Int nPol=lshape(2);
731 0 : posMaxAbs.resize(nPol);
732 0 : for(Int i=0;i<nPol;i++)
733 0 : posMaxAbs(i)=IPosition(lattice.shape().nelements(), 0);
734 0 : maxAbs.resize(nPol);
735 0 : ndx=0;
736 :
737 0 : for(Int s2=0;s2<lshape(2);s2++)
738 0 : for(Int s3=0;s3<lshape(3);s3++)
739 : {
740 0 : ndx(2) = s2; ndx(3)=s3;
741 : {
742 : //
743 : // Locate the pixel with the peak value. That's the
744 : // origin in pixel co-ordinates.
745 : //
746 0 : maxAbs(s2)=0;
747 0 : posMaxAbs(s2) = 0;
748 0 : for(ndx(1)=0;ndx(1)<lshape(1);ndx(1)++)
749 0 : for(ndx(0)=0;ndx(0)<lshape(0);ndx(0)++)
750 0 : if (abs(lattice(ndx)) > maxAbs(s2))
751 0 : {posMaxAbs(s2) = ndx;maxAbs(s2)=abs(lattice(ndx));}
752 : }
753 : }
754 0 : }
755 : //
756 : //---------------------------------------------------------------
757 : // Get the value of the named variable from ~/.aipsrc (or ~/.casarc)
758 : // or from a env. variable (in this precidence order).
759 : //
760 : template <class T>
761 625 : T SynthesisUtils::getenv(const char *name, const T defaultVal)
762 : {
763 625 : T val=defaultVal;
764 625 : stringstream defaultStr;
765 625 : defaultStr << defaultVal;
766 : {
767 625 : char *valStr=NULL;
768 625 : std::string tt(name);
769 : unsigned long pos;
770 1146 : while((pos=tt.find(".")) != tt.npos)
771 521 : tt.replace(pos, 1, "_");
772 :
773 625 : if ((valStr = std::getenv(tt.c_str())) != NULL)
774 : {
775 648 : stringstream toT2(valStr);
776 324 : toT2 >> val;
777 324 : }
778 625 : }
779 : // If environment variable was not found (val remained set to the
780 : // defaulVal), look in ~/.aipsrc.
781 625 : if (val==defaultVal)
782 : {
783 527 : uint handle = Aipsrc::registerRC(name, defaultStr.str().c_str());
784 527 : String strVal = Aipsrc::get(handle);
785 527 : stringstream toT(strVal);
786 527 : toT >> val;
787 527 : }
788 625 : return val;
789 625 : }
790 : template
791 : Int SynthesisUtils::getenv(const char *name, const Int defaultVal);
792 : template
793 : Bool SynthesisUtils::getenv(const char *name, const Bool defaultVal);
794 : template
795 : Float SynthesisUtils::getenv(const char *name, const Float defaultVal);
796 : template
797 : double SynthesisUtils::getenv(const char *name, const double defaultVal);
798 : template
799 : String SynthesisUtils::getenv(const char *name, const String defaultVal);
800 :
801 0 : Float SynthesisUtils::libreSpheroidal(Float nu)
802 : {
803 : Double top, bot, nuend, delnusq;
804 : uInt part;
805 0 : if (nu <= 0) return 1.0;
806 : else
807 0 : if (nu >= 1.0)
808 0 : return 0.0;
809 : else
810 : {
811 0 : uInt np = 5;
812 0 : uInt nq = 3;
813 0 : Matrix<Double> p(np, 2);
814 0 : Matrix<Double> q(nq, 2);
815 0 : p(0,0) = 8.203343e-2;
816 0 : p(1,0) = -3.644705e-1;
817 0 : p(2,0) = 6.278660e-1;
818 0 : p(3,0) = -5.335581e-1;
819 0 : p(4,0) = 2.312756e-1;
820 0 : p(0,1) = 4.028559e-3;
821 0 : p(1,1) = -3.697768e-2;
822 0 : p(2,1) = 1.021332e-1;
823 0 : p(3,1) = -1.201436e-1;
824 0 : p(4,1) = 6.412774e-2;
825 0 : q(0,0) = 1.0000000e0;
826 0 : q(1,0) = 8.212018e-1;
827 0 : q(2,0) = 2.078043e-1;
828 0 : q(0,1) = 1.0000000e0;
829 0 : q(1,1) = 9.599102e-1;
830 0 : q(2,1) = 2.918724e-1;
831 :
832 0 : part = 0;
833 0 : nuend = 0.0;
834 0 : if ((nu >= 0.0) && (nu < 0.75))
835 : {
836 0 : part = 0;
837 0 : nuend = 0.75;
838 : }
839 0 : else if ((nu >= 0.75) && (nu <= 1.00))
840 : {
841 0 : part = 1;
842 0 : nuend = 1.0;
843 : }
844 :
845 0 : top = p(0,part);
846 0 : delnusq = pow(nu,2.0) - pow(nuend,2.0);
847 0 : for (uInt k=1; k<np; k++)
848 0 : top += p(k, part) * pow(delnusq, (Double)k);
849 :
850 0 : bot = q(0, part);
851 0 : for (uInt k=1; k<nq; k++)
852 0 : bot += q(k,part) * pow(delnusq, (Double)k);
853 :
854 0 : if (bot != 0.0) return (top/bot);
855 0 : else return 0.0;
856 0 : }
857 : }
858 0 : Double SynthesisUtils::getRefFreq(const VisBuffer2& /*vb*/)
859 : {
860 0 : throw(AipsError("SynthesisUtils::getRefFreq() depricated due to VI2/VB2 move"));
861 : // return max((vb.getVi()->ms())//vb.msColumns()
862 : // .spectralWindow().chanFreq().getColumn());
863 : }
864 :
865 296 : void SynthesisUtils::makeFTCoordSys(const CoordinateSystem& coords,
866 : const Int& convSize,
867 : const Vector<Double>& ftRef,
868 : CoordinateSystem& ftCoords)
869 : {
870 : Int directionIndex;
871 :
872 296 : ftCoords = coords;
873 296 : directionIndex=ftCoords.findCoordinate(Coordinate::DIRECTION);
874 : // The following line follows the (lame) logic that if a
875 : // DIRECTION axis was not found, the coordinate system must be of
876 : // the FT domain already
877 296 : if (directionIndex == -1) return;
878 :
879 212 : DirectionCoordinate dc;//=coords.directionCoordinate(directionIndex);
880 : // AlwaysAssert(directionIndex>=0, AipsError);
881 212 : dc=coords.directionCoordinate(directionIndex);
882 212 : Vector<Bool> axes(2); axes(0)=axes(1)=true;//axes(2)=true;
883 212 : Vector<Int> shape(2,convSize);
884 :
885 212 : Vector<Double>ref(4);
886 212 : ref(0)=ref(1)=ref(2)=ref(3)=0;
887 212 : dc.setReferencePixel(ref);
888 212 : Coordinate* ftdc=dc.makeFourierCoordinate(axes,shape);
889 212 : Vector<Double> refVal;
890 212 : refVal=ftdc->referenceValue();
891 : // refVal(0)=refVal(1)=0;
892 : // ftdc->setReferenceValue(refVal);
893 212 : ref(0)=ftRef(0);
894 212 : ref(1)=ftRef(1);
895 212 : ftdc->setReferencePixel(ref);
896 :
897 212 : ftCoords.replaceCoordinate(*ftdc, directionIndex);
898 : // LogIO logio;
899 : // IPosition tt;
900 : // coords.list(logio,MDoppler::RADIO,tt,tt);
901 : // ftCoords.list(logio,MDoppler::RADIO,tt,tt);
902 :
903 212 : delete ftdc; ftdc=0;
904 212 : }
905 :
906 : //
907 : // Given a list of Spw,MinFreq,MaxFreq,FreqStep (the output product
908 : // of MSSelection), expand the list to a list of channel freqs. and
909 : // conjugate freqs. per SPW.
910 : //
911 186 : void SynthesisUtils::expandFreqSelection(const Matrix<Double>& freqSelection,
912 : Matrix<Double>& expandedFreqList,
913 : Matrix<Double>& expandedConjFreqList)
914 : {
915 186 : Int nSpw = freqSelection.shape()(0), maxSlots=0;
916 : Double freq;
917 :
918 736 : for (Int s=0;s<nSpw;s++)
919 550 : maxSlots=max(maxSlots,SynthesisUtils::nint((freqSelection(s,2)-freqSelection(s,1))/freqSelection(s,3))+1);
920 :
921 186 : expandedFreqList.resize(nSpw,maxSlots);
922 186 : expandedConjFreqList.resize(nSpw,maxSlots);
923 :
924 736 : for (Int s=0,cs=(nSpw-1);s<nSpw;s++,cs--)
925 1100 : for (Int i=0,ci=(maxSlots-1);i<maxSlots;i++,ci--)
926 : {
927 550 : freq = freqSelection(s,1)+i*freqSelection(s,3);
928 550 : expandedFreqList(s,i) = (freq <= freqSelection(s,2)) ? freq : 0;
929 550 : freq = freqSelection(cs,2) - ci*freqSelection(cs,3);
930 550 : expandedConjFreqList(s,ci) = (freq >= freqSelection(cs,1)) ? freq : 0;
931 : }
932 186 : }
933 :
934 : //
935 : // The result will be in-place in c1
936 : //
937 : template
938 : void SynthesisUtils::libreConvolver(Array<Complex>& c1, const Array<Complex>& c2);
939 :
940 :
941 : template <class T>
942 0 : void SynthesisUtils::libreConvolver(Array<T>& c1, const Array<T>& c2)
943 : {
944 0 : Array<T> c2tmp;
945 0 : c2tmp.assign(c2);
946 :
947 0 : if (c1.shape().product() > c2tmp.shape().product())
948 0 : c2tmp.resize(c1.shape(),true);
949 : else
950 0 : c1.resize(c2tmp.shape(),true);
951 :
952 :
953 0 : ArrayLattice<T> c2tmp_lat(c2tmp), c1_lat(c1);
954 :
955 0 : LatticeFFT::cfft2d(c1_lat,false);
956 0 : LatticeFFT::cfft2d(c2tmp_lat,false);
957 : //cerr << "########## " << c1.shape() << " " << c2tmp.shape() << endl;
958 0 : c1 = sqrt(c1);
959 0 : c2tmp=sqrt(c2tmp);
960 0 : c1 *= conj(c2tmp);
961 0 : LatticeFFT::cfft2d(c1_lat);
962 0 : }
963 :
964 :
965 1160 : Double SynthesisUtils::nearestValue(const Vector<Double>& list, const Double& val, Int& index)
966 : {
967 2320 : Vector<Double> diff = fabs(list - val);
968 1160 : Double minVal=1e20;
969 1160 : Int i=0;
970 :
971 4600 : for (index=0;index<(Int)diff.nelements();index++)
972 3440 : if (diff[index] < minVal) {minVal=diff[index];i=index;}
973 1160 : index=i;
974 2320 : return list(index);
975 :
976 : // The algorithm below has a N*log(N) cost.
977 : //
978 : // Bool dummy;
979 : // Sort sorter(diff.getStorage(dummy), sizeof(Double));
980 : // sorter.sortKey((uInt)0,TpDouble);
981 : // Int nch=list.nelements();
982 : // Vector<uInt> sortIndx;
983 : // sorter.sort(sortIndx, nch);
984 :
985 : // index=sortIndx(0);
986 : // return list(index);
987 :
988 :
989 : // Works only for regularly samplied list
990 : //
991 : // Int ndx=min(freqValues_p.nelements()-1,max(0,SynthesisUtils::nint((freqVal-freqValues_p[0])/freqValIncr_p)));
992 : // return ndx;
993 1160 : }
994 :
995 : template <class T>
996 4464 : T SynthesisUtils::stdNearestValue(const vector<T>& list, const T& val, Int& index)
997 : {
998 : // auto const it = std::lower_bound(list.begin(), list.end(), val);
999 : // if (it == list.begin()) return list[0];
1000 : // else return list[*(it-1)];
1001 :
1002 4464 : vector<T> diff=list;
1003 12276 : for (uInt i=0;i<list.size();i++)
1004 7812 : diff[i] = fabs(list[i] - val);
1005 :
1006 4464 : T minVal=std::numeric_limits<T>::max();//1e20;
1007 4464 : Int i=0;
1008 :
1009 12276 : for (index=0;index<(Int)diff.size();index++)
1010 7812 : if (diff[index] < minVal) {minVal=diff[index];i=index;}
1011 4464 : index=i;
1012 8928 : return list[index];
1013 4464 : }
1014 :
1015 0 : CoordinateSystem SynthesisUtils::makeUVCoords(const CoordinateSystem& imageCoordSys,
1016 : const IPosition& shape)
1017 : {
1018 0 : CoordinateSystem FTCoords = imageCoordSys;
1019 :
1020 0 : Int dirIndex=FTCoords.findCoordinate(Coordinate::DIRECTION);
1021 0 : DirectionCoordinate dc=imageCoordSys.directionCoordinate(dirIndex);
1022 0 : Vector<Bool> axes(2); axes=true;
1023 0 : Vector<Int> dirShape(2); dirShape(0)=shape(0);dirShape(1)=shape(1);
1024 0 : Coordinate* FTdc=dc.makeFourierCoordinate(axes,dirShape);
1025 0 : FTCoords.replaceCoordinate(*FTdc,dirIndex);
1026 0 : delete FTdc;
1027 :
1028 0 : return FTCoords;
1029 0 : }
1030 :
1031 282 : Vector<Int> SynthesisUtils::mapSpwIDToDDID(const VisBuffer2& vb, const Int& spwID)
1032 : {
1033 282 : Vector<Int> ddidList;
1034 : //Int nDDRows = vb.msColumns().dataDescription().nrow();
1035 282 : Int nDDRows = (vb.ms()).dataDescription().nrow();
1036 1124 : for (Int i=0; i<nDDRows; i++)
1037 : {
1038 842 : if((vb.subtableColumns()).dataDescription().spectralWindowId().get(i) == spwID)
1039 : {
1040 282 : Int n=ddidList.nelements();
1041 282 : ddidList.resize(n+1,true);
1042 282 : ddidList(n) = i;
1043 : }
1044 : }
1045 282 : return ddidList;
1046 0 : }
1047 :
1048 282 : Vector<Int> SynthesisUtils::mapSpwIDToPolID(const VisBuffer2& vb, const Int& spwID)
1049 : {
1050 282 : Vector<Int> polIDList, ddIDList;
1051 282 : ddIDList = SynthesisUtils::mapSpwIDToDDID(vb, spwID);
1052 282 : Int n=ddIDList.nelements();
1053 282 : polIDList.resize(n);
1054 564 : for (Int i=0; i<n; i++)
1055 282 : polIDList(i) = (vb.subtableColumns()).dataDescription().polarizationId().get(ddIDList(i));
1056 :
1057 564 : return polIDList;
1058 282 : }
1059 :
1060 :
1061 : //
1062 : // Calcualte the BLC, TRC of the intersection of two rectangles (courtesy U.Rau)
1063 : //
1064 0 : void SynthesisUtils::calcIntersection(const Int blc1[2], const Int trc1[2],
1065 : const Float blc2[2], const Float trc2[2],
1066 : Float blc[2], Float trc[2])
1067 : {
1068 : // cerr << blc1[0] << " " << blc1[1] << " " << trc1[0] << " " << trc1[1] << " " << blc2[0] << " " << blc2[1] << " " << trc2[0] << " " << trc2[1] << endl;
1069 : Float dblc, dtrc;
1070 0 : for (Int i=0;i<2;i++)
1071 : {
1072 0 : dblc = blc2[i] - blc1[i];
1073 0 : dtrc = trc2[i] - trc1[i];
1074 :
1075 0 : if ((dblc >= 0) and (dtrc >= 0))
1076 : {
1077 0 : blc[i] = blc1[i] + dblc;
1078 0 : trc[i] = trc2[i] - dtrc;
1079 : }
1080 0 : else if ((dblc >= 0) and (dtrc < 0))
1081 : {
1082 0 : blc[i] = blc1[i] + dblc;
1083 0 : trc[i] = trc1[i] + dtrc;
1084 : }
1085 0 : else if ((dblc < 0) and (dtrc >= 0))
1086 : {
1087 0 : blc[i] = blc2[i] - dblc;
1088 0 : trc[i] = trc2[i] - dtrc;
1089 : }
1090 : else
1091 : {
1092 0 : blc[i] = blc2[i] - dblc;
1093 0 : trc[i] = trc1[i] + dtrc;
1094 : }
1095 : }
1096 0 : }
1097 : //
1098 : // Check if the two rectangles interset (courtesy U.Rau).
1099 : //
1100 0 : Bool SynthesisUtils::checkIntersection(const Int blc1[2], const Int trc1[2], const Float blc2[2], const Float trc2[2])
1101 : {
1102 : // blc1[2] = {xmin1, ymin1};
1103 : // blc2[2] = {xmin2, ymin2};
1104 : // trc1[2] = {xmax1, ymax1};
1105 : // trc2[2] = {xmax2, ymax2};
1106 :
1107 0 : if ((blc1[0] > trc2[0]) || (trc1[0] < blc2[0]) || (blc1[1] > trc2[1]) || (trc1[1] < blc2[1]))
1108 0 : return false;
1109 : else
1110 0 : return true;
1111 : // def checkintersect( xmin1, ymin1, xmax1, ymax1, xmin2, ymin2, xmax2, ymax2 ):
1112 : // if xmin1 > xmax2 or xmax1 < xmin2 or ymin1 > ymax2 or ymax1 < ymin2 :
1113 : // return false
1114 : // else :
1115 : // return true
1116 : }
1117 :
1118 : // template<class Iterator>
1119 : // Iterator SynthesisUtils::Unique(Iterator first, Iterator last)
1120 : // {
1121 : // while (first != last)
1122 : // {
1123 : // Iterator next(first);
1124 : // last = std::remove(++next, last, *first);
1125 : // first = next;
1126 : // }
1127 :
1128 : // return last;
1129 : // }
1130 :
1131 0 : String SynthesisUtils::mjdToString(casacore::Time& time)
1132 : {
1133 0 : String tStr;
1134 0 : tStr = String::toString(time.year()) + "/" +
1135 0 : String::toString(time.month()) + "/" +
1136 0 : String::toString(time.dayOfMonth()) + "/" +
1137 0 : String::toString(time.hours()) + ":" +
1138 0 : String::toString(time.minutes()) + ":";
1139 0 : ostringstream fsec;
1140 0 : fsec << setprecision(2) << time.dseconds();
1141 0 : tStr = tStr + String(fsec.str());
1142 : // String::toString(time.dseconds());
1143 0 : return tStr;
1144 0 : }
1145 :
1146 : template <class Iterator>
1147 558 : Iterator SynthesisUtils::Unique(Iterator first, Iterator last)
1148 : {
1149 1674 : while (first != last)
1150 : {
1151 1116 : Iterator next(first);
1152 1116 : last = std::remove(++next, last, *first);
1153 1116 : first = next;
1154 : }
1155 :
1156 558 : return last;
1157 : }
1158 :
1159 0 : void SynthesisUtils::showCS(const CoordinateSystem& cs, std::ostream &os, const String& msg)
1160 : {
1161 0 : LogIO log_l(LogOrigin("SynthesisUtils","showCS"));
1162 0 : IPosition dummy;
1163 0 : Vector<String> csList;
1164 0 : if (msg!="")
1165 0 : os << "CoordSys: ";
1166 0 : csList = cs.list(log_l,MDoppler::RADIO,dummy,dummy);
1167 0 : os << csList << std::endl;
1168 0 : }
1169 968 : const casacore::Array<Complex> SynthesisUtils::getCFPixels(const casacore::String& Dir,
1170 : const casacore::String& fileName)
1171 : {
1172 : try
1173 : {
1174 1936 : casacore::PagedImage<casacore::Complex> thisCF(Dir+'/'+fileName);
1175 968 : return thisCF.get();
1176 968 : }
1177 0 : catch (AipsError &x)
1178 : {
1179 0 : LogIO log_l(LogOrigin("SynthesisUtils","getCFPixels"));
1180 0 : log_l << x.getMesg() << LogIO::EXCEPTION;
1181 0 : }
1182 0 : return casacore::Array<Complex>(); // Just to keep the complier happy. Program control should never get here.
1183 : }
1184 :
1185 0 : void SynthesisUtils::putCFPixels(const casacore::String& Dir,
1186 : const casacore::String& fileName,
1187 : const casacore::Array<Complex>& srcpix)
1188 : {
1189 : try
1190 : {
1191 0 : casacore::PagedImage<casacore::Complex> thisCF(Dir+'/'+fileName);
1192 0 : return thisCF.put(srcpix);
1193 0 : }
1194 0 : catch (AipsError &x)
1195 : {
1196 0 : LogIO log_l(LogOrigin("SynthesisUtils","putCFPixels"));
1197 0 : log_l << x.getMesg() << LogIO::EXCEPTION;
1198 0 : }
1199 : }
1200 :
1201 716 : const casacore::IPosition SynthesisUtils::getCFShape(const casacore::String& Dir,
1202 : const casacore::String& fileName)
1203 : {
1204 : try
1205 : {
1206 1432 : casacore::PagedImage<casacore::Complex> thisCF(Dir+'/'+fileName);
1207 716 : return thisCF.shape();
1208 716 : }
1209 0 : catch (AipsError &x)
1210 : {
1211 0 : LogIO log_l(LogOrigin("SynthesisUtils","getCFShape"));
1212 0 : log_l << x.getMesg() << LogIO::EXCEPTION;
1213 0 : }
1214 0 : return casacore::IPosition(); // Just to keep the complier happy. Program control should never get here.
1215 : }
1216 :
1217 2232 : casacore::TableRecord SynthesisUtils::getCFParams(const casacore::String& Dir,
1218 : const casacore::String& fileName,
1219 : casacore::IPosition& cfShape,
1220 : casacore::Array<Complex>& pixelBuffer,
1221 : casacore::CoordinateSystem& coordSys,
1222 : casacore::Double& sampling,
1223 : casacore::Double& paVal,
1224 : casacore::Int& xSupport, casacore::Int& ySupport,
1225 : casacore::Double& fVal, casacore::Double& wVal, casacore::Int& mVal,
1226 : casacore::Double& conjFreq, casacore::Int& conjPoln,
1227 : casacore::Bool loadPixels,
1228 : casacore::Bool loadMiscInfo)
1229 : {
1230 : try
1231 : {
1232 : //casacore::Table tThisCF(Dir+'/'+fileName);
1233 4464 : casacore::PagedImage<casacore::Complex> thisCF(Dir+'/'+fileName);
1234 2232 : cfShape = thisCF.shape();
1235 2232 : if (loadPixels) pixelBuffer.assign(thisCF.get());
1236 2232 : casacore::TableRecord miscinfo;
1237 2232 : if (loadMiscInfo)
1238 : {
1239 2232 : miscinfo= thisCF.miscInfo();
1240 :
1241 2232 : miscinfo.get("ParallacticAngle", paVal);
1242 2232 : miscinfo.get("MuellerElement", mVal);
1243 2232 : miscinfo.get("WValue", wVal);
1244 2232 : miscinfo.get("Xsupport", xSupport);
1245 2232 : miscinfo.get("Ysupport", ySupport);
1246 2232 : miscinfo.get("Sampling", sampling);
1247 2232 : miscinfo.get("ConjFreq", conjFreq);
1248 2232 : miscinfo.get("ConjPoln", conjPoln);
1249 2232 : casacore::Int index= thisCF.coordinates().findCoordinate(casacore::Coordinate::SPECTRAL);
1250 2232 : coordSys = thisCF.coordinates();
1251 2232 : casacore::SpectralCoordinate spCS = coordSys.spectralCoordinate(index);
1252 2232 : fVal=static_cast<Float>(spCS.referenceValue()(0));
1253 2232 : }
1254 4464 : return miscinfo;
1255 2232 : }
1256 0 : catch(AipsError& x)
1257 : {
1258 0 : throw(AipsError(String("Error in SynthesisUtils::getCFParams(): ")
1259 0 : +x.getMesg()));
1260 0 : }
1261 : };
1262 :
1263 :
1264 920 : void SynthesisUtils::rotate2(const double& actualPA, CFCell& baseCFC, CFCell& cfc, const double& rotAngleIncr)
1265 : {
1266 1840 : LogIO log_l(LogOrigin("SynthesisUtils", "rotate2"));
1267 :
1268 : // // If the A-Term is a No-Op, the resulting CF is rotationally
1269 : // // symmetric.
1270 : // if (isNoOp()) return;
1271 :
1272 : (void)baseCFC;
1273 :
1274 : //double actualPA = getPA(vb);
1275 920 : double currentCFPA = cfc.pa_p.getValue("rad");
1276 : //Double baseCFCPA=baseCFC.pa_p.getValue("rad");
1277 :
1278 920 : double dPA = currentCFPA-actualPA;
1279 :
1280 920 : if (fabs(dPA) > fabs(rotAngleIncr))
1281 : {
1282 0 : casacore::Array<TT> inData;
1283 : //inData.assign(*baseCFC.getStorage());
1284 : //dPA = baseCFCPA-actualPA;
1285 0 : dPA = currentCFPA-actualPA;
1286 0 : inData.assign(*cfc.getStorage());
1287 : try
1288 : {
1289 0 : SynthesisUtils::rotateComplexArray(log_l, inData, cfc.coordSys_p,
1290 0 : *cfc.getStorage(),
1291 : dPA);//,"LINEAR");
1292 : // currentCFPA-actualPA);//,"LINEAR");
1293 : }
1294 0 : catch (AipsError &x)
1295 : {
1296 0 : log_l << x.getMesg() << LogIO::EXCEPTION;
1297 0 : }
1298 0 : cfc.pa_p=Quantity(actualPA, "rad");
1299 0 : }
1300 920 : };
1301 :
1302 :
1303 : //
1304 : // Parser for parsing the NAME field of the SPECTRAL_WINDOW sub-table.
1305 : // Returns a vector of string containing tokens separated by "#"
1306 : // in the supplied band name string.
1307 : //
1308 3252 : casacore::Vector<casacore::String> SynthesisUtils::parseBandName(const casacore::String& bandName)
1309 : {
1310 3252 : casacore::Vector<casacore::String> tokens;
1311 :
1312 : // Allow a blank band name for older MSes where this is sometimes left blank.
1313 3252 : if (bandName == "")
1314 : {
1315 1092 : tokens.resize(1);tokens="";
1316 1092 : return tokens;
1317 : }
1318 :
1319 : char *tok, *name;
1320 2160 : int i=0;
1321 :
1322 2160 : name = (char *)bandName.c_str();
1323 2160 : if ((tok = std::strtok(name,"#"))!=NULL)
1324 : {
1325 2160 : tokens.resize(i+1,true); tokens(i)="";
1326 2160 : tokens(i++).assign(tok);
1327 : }
1328 :
1329 2160 : while ((tok = std::strtok(NULL,"#"))!=NULL)
1330 : {
1331 0 : tokens.resize(i+1,true); tokens(i)="";
1332 0 : tokens(i++).assign(tok);
1333 : }
1334 2160 : if (tokens(0)=="")
1335 : {
1336 0 : LogIO log_l(LogOrigin("SynthesisUtils", "parseBandName"));
1337 0 : log_l << "Error while parsing band name \"" << bandName << "\"" << LogIO::EXCEPTION;
1338 0 : }
1339 :
1340 : // for (i=0;i<3;i++)
1341 : // if (tokens(i)=="") log_l << "Error while parsing band name \"" << bandName << LogIO::EXCEPTION;
1342 2160 : return tokens;
1343 0 : }
1344 :
1345 :
1346 : //
1347 : //-----------------------------------------------------------------------------------------
1348 : //
1349 0 : casacore::CoordinateSystem SynthesisUtils::makeModelGridFromImage(const std::string& modelImageName,
1350 : casacore::TempImage<casacore::DComplex>& modelImageGrid)
1351 : {
1352 : // This code is basically loading a floating point image from
1353 : // the disk, and loading it into a complex<double> image. This
1354 : // should be possible on-the-fly.
1355 : //
1356 : // However currently it is not possible to this OTF. So one has
1357 : // to load the image from the disk in a float image (copy-1 of
1358 : // the image in the memory). Then, since
1359 : // StokesImageUtil::From() does not work with complex<double>
1360 : // images, convert it first to a complex<float> image (equal to
1361 : // 2 more float buffers in the memory). And then covert the
1362 : // complex<float> image to a complex<Double> image (equal to 4
1363 : // more float buffers in the memory).
1364 : //
1365 : // In the end, just because of limitations of OTF type
1366 : // conversions, we end up with 7x memory footprint!
1367 :
1368 0 : casacore::LatticeBase* lattPtr = casacore::ImageOpener::openImage (modelImageName);
1369 : casacore::ImageInterface<float> *fImage;
1370 0 : fImage = dynamic_cast<casacore::ImageInterface<float>*>(lattPtr);
1371 :
1372 0 : TempImage<casacore::Complex> tmp(fImage->shape(), fImage->coordinates());
1373 0 : StokesImageUtil::From(tmp, *fImage);
1374 :
1375 0 : modelImageGrid = casacore::TempImage<casacore::DComplex> (fImage->shape(), fImage->coordinates());
1376 :
1377 : Bool d0,d1;
1378 0 : casacore::Array<casacore::DComplex> dcArray=modelImageGrid.get();
1379 0 : casacore::Array<casacore::Complex> fcArray=tmp.get();
1380 :
1381 0 : casacore::DComplex* dcArrayPtr= dcArray.getStorage(d0);
1382 0 : casacore::Complex* fcArrayPtr = fcArray.getStorage(d1);
1383 0 : IPosition ndx(4,0,0,0,0),shape=fImage->shape();
1384 :
1385 0 : for (ndx(0)=0; ndx(0)<shape(0);ndx(0)++)
1386 0 : for (ndx(1)=0; ndx(1)<shape(1);ndx(1)++)
1387 0 : for (ndx(2)=0; ndx(2)<shape(2);ndx(2)++)
1388 0 : for (ndx(3)=0; ndx(3)<shape(3);ndx(3)++)
1389 0 : dcArray(ndx) = DComplex(fcArray(ndx).real(), fcArray(ndx).imag());
1390 :
1391 0 : modelImageGrid.put(dcArray);
1392 0 : return fImage->coordinates();
1393 0 : }
1394 : //
1395 : //-------------------------------------------------------------------------------------
1396 : //
1397 0 : void SynthesisUtils::makeAWLists(const casacore::Vector<double>& wVals,
1398 : const casacore::Vector<double>& fVals,
1399 : const bool& wbAWP, const uint& nw,
1400 : const double& imRefFreq, const double& spwRefFreq,
1401 : casacore::Vector<int>& wNdxList,
1402 : casacore::Vector<int>& spwNdxList,
1403 : const int vbSPW)
1404 : {
1405 : //
1406 : // The following can be generalized to pick a subset of CFs along
1407 : // W and SPW axis in the CFB. Perhaps also useful in the longer
1408 : // run, e.g. when using a super-set CFC not all of which may be
1409 : // used for a given imaging.
1410 : //
1411 : // W-pixels in the CFC should be >= w-planes user setting
1412 0 : assert(wVals.nelements() >= nw);
1413 :
1414 : // Make list of W-CF indexes
1415 0 : int nWCFs=(nw<=1)?nw:wVals.nelements();
1416 0 : wNdxList.resize(nWCFs);
1417 0 : for(int i=0;i<nWCFs;i++) wNdxList[i] = i;
1418 :
1419 : // Make list of SPW-CF indexes
1420 0 : int nSPWCFs=fVals.nelements();
1421 0 : if (wbAWP==true)
1422 : {
1423 : // If a valid SPW ID is given, translate it to the spwNdx for the nearest SPW
1424 0 : if ((vbSPW>=0))// && (vbSPW <nSPWCFs))
1425 : {
1426 : int refSPW;
1427 0 : std::vector<double> stdList(nSPWCFs);
1428 0 : for (int i=0; i<nSPWCFs; i++) stdList[i] = fVals[i];
1429 : //Double refCFFreq =
1430 0 : SynthesisUtils::stdNearestValue(stdList, spwRefFreq, refSPW);
1431 :
1432 0 : spwNdxList.resize(1);
1433 0 : spwNdxList[0]=refSPW;
1434 0 : }
1435 : else
1436 : {
1437 0 : spwNdxList.resize(nSPWCFs);
1438 0 : for(int i=0;i<nSPWCFs;i++) spwNdxList[i] = i;
1439 : }
1440 : }
1441 : else
1442 : {
1443 : // For wbAWP=F, pick up the CF closest to the image reference frequency
1444 : int refSPW;
1445 0 : std::vector<double> stdList(nSPWCFs);
1446 0 : for (int i=0; i<nSPWCFs; i++) stdList[i] = fVals[i];
1447 : //Double refCFFreq =
1448 0 : SynthesisUtils::stdNearestValue(stdList, imRefFreq, refSPW);
1449 :
1450 0 : spwNdxList.resize(1);
1451 0 : spwNdxList[0]=refSPW;
1452 0 : }
1453 :
1454 0 : return;
1455 : }
1456 :
1457 :
1458 : template
1459 : std::vector<Double>::iterator SynthesisUtils::Unique(std::vector<Double>::iterator first, std::vector<Double>::iterator last);
1460 : template
1461 : std::vector<Float>::iterator SynthesisUtils::Unique(std::vector<Float>::iterator first, std::vector<Float>::iterator last);
1462 : template
1463 : std::vector<Int>::iterator SynthesisUtils::Unique(std::vector<Int>::iterator first, std::vector<Int>::iterator last);
1464 :
1465 : template
1466 : Double SynthesisUtils::stdNearestValue(const vector<Double>& list, const Double& val, Int& index);
1467 : template
1468 : Float SynthesisUtils::stdNearestValue(const vector<Float>& list, const Float& val, Int& index);
1469 : template
1470 : Int SynthesisUtils::stdNearestValue(const vector<Int>& list, const Int& val, Int& index);
1471 :
1472 : /////===========
1473 135 : MathUtils::MathUtils(){
1474 : //Timer tim;
1475 : //tim.mark();
1476 135 : initSincCache();
1477 : //tim.show("Calculating 16000 sines");
1478 135 : }
1479 0 : Array<Complex> MathUtils::resample(const Array<Complex>& inarray, const Double factorX, const Double factorY) {
1480 :
1481 0 : if(factorX == 1.0 && factorY==1.0)
1482 0 : return inarray;
1483 0 : Double nx=Double(inarray.shape()(0));
1484 0 : Double ny=Double(inarray.shape()(1));
1485 0 : IPosition shp=inarray.shape();
1486 : //cerr << "shp " << shp << endl;
1487 0 : shp(0)=Int(nx*factorX/8.0)*8;
1488 0 : shp(1)=Int(ny*factorY/8.0)*8;
1489 0 : Int newNx=shp(0);
1490 0 : Int newNy=shp(1);
1491 : // cerr << "SHP " << shp << endl;
1492 0 : Array<Complex> out(shp, Complex(0.0));
1493 :
1494 : /*IPosition incursor=IPosition(inarray.shape().nelements(),1);
1495 : incursor[0]=nx;
1496 : incursor[1]=ny;
1497 : IPosition outcursor=IPosition(inarray.shape().nelements(),1);
1498 : outcursor[0]=shp[0];
1499 : outcursor[1]=shp[1];
1500 : */
1501 0 : ArrayIterator<Complex> inIt(inarray, IPosition(2,0,1), True);
1502 0 : ArrayIterator<Complex> outIt(out, IPosition(2,0,1),True);
1503 0 : inIt.origin();
1504 0 : outIt.origin();
1505 : //for (zzz=0; zzz< shp.(4); ++zzz){
1506 : // for(yyy=0; yyy< shp.(3); ++yyy){
1507 : // for(xxx=0; xxx< shp.(2); ++xxx){
1508 0 : while(!inIt.pastEnd()) {
1509 : // cerr << "Iter shape " << inIt.array().shape() << endl;
1510 0 : Matrix<Complex> inmat;
1511 0 : inmat=inIt.array();
1512 : //Matrix<Float> leReal=real(Matrix<Complex>(inIt.array()));
1513 : //Matrix<Float> leImag=imag(Matrix<Complex>(inIt.array()));
1514 0 : Matrix<Float> leReal=real(inmat);
1515 0 : Matrix<Float> leImag=imag(inmat);
1516 : Bool leRealCopy, leImagCopy;
1517 0 : Float *realptr=leReal.getStorage(leRealCopy);
1518 0 : Float *imagptr=leImag.getStorage(leImagCopy);
1519 : Bool isCopy;
1520 0 : Matrix<Complex> outMat(outIt.array());
1521 0 : Complex *intPtr=outMat.getStorage(isCopy);
1522 : Float realval, imagval;
1523 : #ifdef _OPENMP
1524 0 : omp_set_nested(0);
1525 : #endif
1526 0 : #pragma omp parallel for default(none) private(realval, imagval) firstprivate(intPtr, realptr, imagptr, nx, ny, newNx, newNy) shared(leReal, leImag)
1527 :
1528 : for (Int k =0; k < newNy; ++k) {
1529 : Double y =Double(k)/Double(newNy)*Double(ny);
1530 :
1531 : for (Int j=0; j < newNx; ++j) {
1532 : // Interpolate2D interp(Interpolate2D::LANCZOS);
1533 : Double x=Double(j)/Double(newNx)*Double(nx);
1534 : //interp.interp(realval, where, leReal);
1535 : realval=interpLanczos(x , y, nx, ny,
1536 : realptr, 3);
1537 : imagval=interpLanczos(x , y, nx, ny,
1538 : imagptr, 3);
1539 : //interp.interp(imagval, where, leImag);
1540 : intPtr[k*Int(newNx)+j]=Complex(realval, imagval);
1541 : }
1542 :
1543 : }
1544 0 : outMat.putStorage(intPtr, isCopy);
1545 0 : leReal.putStorage(realptr, leRealCopy);
1546 0 : leImag.putStorage(imagptr, leImagCopy);
1547 0 : inIt.next();
1548 0 : outIt.next();
1549 0 : }
1550 0 : return out;
1551 0 : }
1552 135 : void MathUtils::initSincCache(){
1553 1080135 : for (Float u=-4000; u<4000; ++u){
1554 1080000 : Float ux=u/1000.0;
1555 1080000 : if (ux == 0) {
1556 135 : sincCache_p[u+4000]=1.0;
1557 : }
1558 : else{
1559 1079865 : sincCache_p[u+4000]= sin(C::pi * ux) / (C::pi * ux);
1560 : }
1561 : }
1562 135 : sincCachePtr_p=sincCache_p.data();
1563 :
1564 135 : }
1565 0 : Float MathUtils::sinc(const Float x) {
1566 0 : Int index=x*1000+4000;
1567 :
1568 :
1569 0 : return sincCachePtr_p[index];
1570 :
1571 : }
1572 0 : casacore::Float MathUtils::interpLanczos( const casacore::Double& x , const casacore::Double& y, const casacore::Double& nx, const casacore::Double& ny, const casacore::Float* data, const casacore::Float a){
1573 0 : Double floorx = floor(x);
1574 0 : Double floory = floor(y);
1575 0 : Float result=0.0;
1576 0 : if (floorx < a || floorx >= nx - a || floory < a || floory >= ny - a) {
1577 0 : result = 0;
1578 0 : return result;
1579 : }
1580 0 : for (Float i = floorx - a + 1; i <= floorx + a; ++i) {
1581 0 : for (Float j = floory - a + 1; j <= floory + a; ++j) {
1582 0 : result += Float(Double(data[Int(j*nx+i)]) * sinc(x - i)*sinc((x-i)/ a) * sinc(y - j)*sinc((y-j)/ a));
1583 : }
1584 : }
1585 0 : return result;
1586 :
1587 :
1588 : }
1589 0 : Array<Complex> MathUtils::getMiddle(const Array<Complex>& inArr, const int nx, const int ny){
1590 0 : IPosition outshape=inArr.shape();
1591 0 : outshape[0]=nx;
1592 0 : outshape[1]=ny;
1593 0 : IPosition blc(2, (inArr.shape()[0]-nx)/2, (inArr.shape()[1]-ny)/2);
1594 0 : IPosition trc(2, (inArr.shape()[0]+nx)/2-1, (inArr.shape()[1]+ny)/2-1);
1595 0 : Array<Complex> outArr(outshape);
1596 0 : ArrayIterator<Complex> inIt(inArr, IPosition(2,0,1));
1597 0 : ArrayIterator<Complex> outIt(outArr, IPosition(2,0,1));
1598 0 : inIt.origin();
1599 0 : outIt.origin();
1600 0 : while(!inIt.pastEnd()){
1601 : //cerr << "Shapes in getM " << outIt.array().shape() << " in " << inIt.array()(blc, trc).shape() << endl;
1602 0 : outIt.array().assign(inIt.array()(blc, trc));
1603 0 : inIt.next();
1604 0 : outIt.next();
1605 :
1606 :
1607 : }
1608 :
1609 0 : return outArr;
1610 0 : }
1611 :
1612 0 : void MathUtils::putMiddle(Array<Complex>& outArr, const Array<Complex>& inArr) {
1613 0 : Int nx = inArr.shape()[0];
1614 0 : Int ny = inArr.shape()[1];
1615 0 : if(nx < outArr.shape()[0] && ny < outArr.shape()[1]){
1616 0 : IPosition blc(2, (outArr.shape()[0]-nx)/2, (outArr.shape()[1]-ny)/2);
1617 0 : IPosition trc(2, (outArr.shape()[0]+nx)/2-1, (outArr.shape()[1]+ny)/2-1);
1618 0 : ArrayIterator<Complex> inIt(inArr, IPosition(2,0,1));
1619 0 : ArrayIterator<Complex> outIt(outArr, IPosition(2,0,1));
1620 0 : inIt.origin();
1621 0 : outIt.origin();
1622 0 : while(!inIt.pastEnd() && !inIt.pastEnd()){
1623 :
1624 0 : (outIt.array())(blc, trc).assign(inIt.array());
1625 :
1626 0 : inIt.next();
1627 0 : outIt.next();
1628 : }
1629 0 : }
1630 0 : else if(outArr.shape()[0] < nx && outArr.shape()[1] < ny){// take the inner of inArray
1631 0 : IPosition blc(2, (nx-outArr.shape()[0])/2, (ny-outArr.shape()[1])/2);
1632 0 : IPosition trc(2, (outArr.shape()[0]+nx)/2-1, (outArr.shape()[1]+ny)/2-1);
1633 0 : ArrayIterator<Complex> inIt(inArr, IPosition(2,0,1));
1634 0 : ArrayIterator<Complex> outIt(outArr, IPosition(2,0,1));
1635 0 : inIt.origin();
1636 0 : outIt.origin();
1637 0 : while(!inIt.pastEnd() && !inIt.pastEnd()){
1638 : //cerr << "Shapes in putM " << outIt.array().shape() << " in " << inIt.array()(blc, trc).shape() << endl;
1639 0 : outIt.array()=inIt.array()(blc, trc);
1640 0 : inIt.next();
1641 0 : outIt.next();
1642 : }
1643 :
1644 :
1645 0 : }
1646 : else{
1647 0 : throw(AipsError("Programmer's error cannot use PutMiddle"));
1648 :
1649 : }
1650 :
1651 :
1652 0 : }
1653 :
1654 0 : Array<Complex> MathUtils::resampleViaFFT(const Array<Complex>& inarray, const Double factorX, const Double factorY) {
1655 :
1656 0 : if(factorX==1.0 && factorY==1.0)
1657 0 : return inarray;
1658 0 : Double nx=Double(inarray.shape()(0));
1659 0 : Double ny=Double(inarray.shape()(1));
1660 0 : IPosition shp=inarray.shape();
1661 : //cerr << "shp " << shp << endl;
1662 0 : shp(0)=Int(std::ceil(nx*factorX/8.0))*8;
1663 0 : shp(1)=Int(std::ceil(ny*factorY/8.0))*8;
1664 0 : Int newNx=shp(0);
1665 0 : Int newNy=shp(1);
1666 : /* cerr << "SHP " << shp << endl;
1667 : Array<Complex> out(shp, Complex(0.0));
1668 : ArrayIterator<Complex> inIt(inarray, IPosition(2,0,1), True);
1669 : ArrayIterator<Complex> outIt(out, IPosition(2,0,1),True);
1670 : inIt.origin();
1671 : outIt.origin();
1672 : FFT2D ftsmall;
1673 : FFT2D ftlarge;
1674 :
1675 : while(!inIt.pastEnd()) {
1676 : // cerr << "Iter shape " << inIt.array().shape() << endl;
1677 : Matrix<Complex> inmat;
1678 : inmat=inIt.array();
1679 : Bool isCopy;
1680 : Complex * inmatptr=inmat.getStorage(isCopy);
1681 : ftsmall.c2cFFT(inmatptr, nx, ny, True);
1682 : inmat.putStorage(inmatptr,isCopy);
1683 : Matrix<Complex> outMat(outIt.array());
1684 : putMiddle(outMat, inmat);
1685 : Complex *intPtr=outMat.getStorage(isCopy);
1686 : ftlarge.c2cFFT(intPtr, newNx, newNy, False);
1687 : outMat.putStorage(intPtr, isCopy);
1688 : Float fac=Float(newNx)*Float(newNy)/Float(nx*ny);
1689 : outMat *= fac;
1690 : inIt.next();
1691 : outIt.next();
1692 :
1693 : }*/
1694 :
1695 0 : return resampleViaFFT(inarray, newNx, newNy);
1696 0 : }
1697 0 : Array<Complex> MathUtils::resampleViaFFT(const Array<Complex>& inarray, const Int newNx, const Int newNy) {
1698 :
1699 :
1700 0 : Double nx=Double(inarray.shape()(0));
1701 0 : Double ny=Double(inarray.shape()(1));
1702 0 : if (newNx == nx && newNy == ny)
1703 0 : return inarray;
1704 0 : IPosition shp=inarray.shape();
1705 0 : cerr << "shp " << shp << endl;
1706 :
1707 0 : shp(0) = newNx;
1708 0 : shp(1) = newNy;
1709 0 : cerr << "SHP " << shp << endl;
1710 0 : Array<Complex> out(shp, Complex(0.0));
1711 0 : ArrayIterator<Complex> inIt(inarray, IPosition(2,0,1), True);
1712 0 : ArrayIterator<Complex> outIt(out, IPosition(2,0,1),True);
1713 0 : inIt.origin();
1714 0 : outIt.origin();
1715 0 : FFT2D ftsmall;
1716 0 : FFT2D ftlarge;
1717 :
1718 0 : while(!inIt.pastEnd()) {
1719 : // cerr << "Iter shape " << inIt.array().shape() << endl;
1720 0 : Matrix<Complex> inmat;
1721 0 : inmat=inIt.array();
1722 : Bool isCopy;
1723 0 : Complex * inmatptr=inmat.getStorage(isCopy);
1724 0 : ftsmall.c2cFFT(inmatptr, nx, ny, True);
1725 0 : inmat.putStorage(inmatptr,isCopy);
1726 0 : Matrix<Complex> outMat(outIt.array());
1727 0 : putMiddle(outMat, inmat);
1728 0 : Complex *intPtr=outMat.getStorage(isCopy);
1729 0 : ftlarge.c2cFFT(intPtr, newNx, newNy, False);
1730 0 : outMat.putStorage(intPtr, isCopy);
1731 0 : Float fac=Float(newNx)*Float(newNy)/Float(nx*ny);
1732 0 : outMat *= fac;
1733 0 : inIt.next();
1734 0 : outIt.next();
1735 :
1736 0 : }
1737 :
1738 0 : return out;
1739 0 : }
1740 : }
1741 :
1742 : //using namespace casacore;
1743 : } // namespace casa
|
