Line data Source code
1 : //# SetJyGridFT.cc: Implementation of GridFT class
2 : //# Copyright (C) 2012
3 : //# Associated Universities, Inc. Washington DC, USA.
4 : //#
5 : //# This library is free software; you can redistribute it and/or modify it
6 : //# under the terms of the GNU Library General Public License as published by
7 : //# the Free Software Foundation; either version 2 of the License, or (at your
8 : //# option) any later version.
9 : //#
10 : //# This library is distributed in the hope that it will be useful, but WITHOUT
11 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
13 : //# License for more details.
14 : //#
15 : //# You should have received a copy of the GNU Library General Public License
16 : //# along with this library; if not, write to the Free Software Foundation,
17 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
18 : //#
19 : //# Correspondence concerning AIPS++ should be addressed as follows:
20 : //# Internet email: aips2-request@nrao.edu.
21 : //# Postal address: AIPS++ Project Office
22 : //# National Radio Astronomy Observatory
23 : //# 520 Edgemont Road
24 : //# Charlottesville, VA 22903-2475 USA
25 : //#
26 : /*
27 : * SetJyGridFT.cc
28 : *
29 : * Created on: Jun 11, 2012
30 : * Author: kgolap
31 : */
32 : #include <msvis/MSVis/VisibilityIterator.h>
33 : #include <casacore/casa/Quanta/UnitMap.h>
34 : #include <casacore/casa/Quanta/UnitVal.h>
35 : #include <casacore/measures/Measures/Stokes.h>
36 : #include <casacore/coordinates/Coordinates/CoordinateSystem.h>
37 : #include <casacore/coordinates/Coordinates/DirectionCoordinate.h>
38 : #include <casacore/coordinates/Coordinates/SpectralCoordinate.h>
39 : #include <casacore/coordinates/Coordinates/StokesCoordinate.h>
40 : #include <casacore/coordinates/Coordinates/Projection.h>
41 : #include <casacore/ms/MeasurementSets/MSColumns.h>
42 : #include <casacore/casa/BasicSL/Constants.h>
43 : #include <casacore/scimath/Mathematics/FFTServer.h>
44 : #include <synthesis/TransformMachines/SetJyGridFT.h>
45 : #include <casacore/scimath/Mathematics/RigidVector.h>
46 : #include <msvis/MSVis/StokesVector.h>
47 : #include <synthesis/TransformMachines/StokesImageUtil.h>
48 : #include <msvis/MSVis/VisBuffer.h>
49 : #include <msvis/MSVis/VisSet.h>
50 : #include <casacore/images/Images/ImageInterface.h>
51 : #include <casacore/images/Images/PagedImage.h>
52 : #include <casacore/casa/Containers/Block.h>
53 : #include <casacore/casa/Containers/Record.h>
54 : #include <casacore/casa/Arrays/ArrayLogical.h>
55 : #include <casacore/casa/Arrays/ArrayMath.h>
56 : #include <casacore/casa/Arrays/Array.h>
57 : #include <casacore/casa/Arrays/MaskedArray.h>
58 : #include <casacore/casa/Arrays/Vector.h>
59 : #include <casacore/casa/Arrays/Matrix.h>
60 : #include <casacore/casa/Arrays/Cube.h>
61 : #include <casacore/casa/Arrays/MatrixIter.h>
62 : #include <casacore/casa/BasicSL/String.h>
63 : #include <casacore/casa/Utilities/Assert.h>
64 : #include <casacore/casa/Exceptions/Error.h>
65 : #include <casacore/lattices/Lattices/ArrayLattice.h>
66 : #include <casacore/measures/Measures/UVWMachine.h>
67 : #include <casacore/lattices/Lattices/SubLattice.h>
68 : #include <casacore/lattices/LRegions/LCBox.h>
69 : #include <casacore/lattices/Lattices/LatticeCache.h>
70 : #include <casacore/lattices/LatticeMath/LatticeFFT.h>
71 : #include <casacore/lattices/Lattices/LatticeIterator.h>
72 : #include <casacore/lattices/Lattices/LatticeStepper.h>
73 : #include <casacore/scimath/Mathematics/ConvolveGridder.h>
74 : #include <casacore/casa/Utilities/CompositeNumber.h>
75 : #include <casacore/casa/OS/Timer.h>
76 : #include <sstream>
77 : #ifdef _OPENMP
78 : #include <omp.h>
79 : #endif
80 :
81 : using namespace casacore;
82 : namespace casa { //# NAMESPACE CASA - BEGIN
83 : using namespace casacore;
84 : // using namespace casa::async;
85 :
86 : #define NEED_UNDERSCORES
87 : #if defined(NEED_UNDERSCORES)
88 : #define sectdgridjy sectdgridjy_
89 : #endif
90 :
91 : extern "C" {
92 : void sectdgridjy(Complex*,
93 : const Int*,
94 : const Int*,
95 : const Int*,
96 : const Int*,
97 : const Int*,
98 : const Complex*,
99 : const Int*,
100 : const Int*,
101 : const Int *,
102 : const Int *,
103 : //support
104 : const Int*,
105 : const Int*,
106 : const Double*,
107 : const Int*,
108 : const Int*,
109 : //rbeg, rend, loc, off, phasor
110 : const Int*,
111 : const Int*,
112 : const Int*,
113 : const Int*,
114 : const Complex*,
115 : const Double*);
116 : }
117 4 : SetJyGridFT::SetJyGridFT(Long icachesize, Int itilesize, String iconvType,
118 : MPosition mLocation, MDirection mTangent, Float padding,
119 4 : Bool usezero, Bool useDoublePrec,const Vector<Double>& frequencyscale, const Vector<Double>& scale):GridFT(icachesize, itilesize, iconvType,
120 : mLocation, mTangent, padding, usezero, useDoublePrec), freqscale_p(frequencyscale),
121 4 : scale_p(scale) {
122 :
123 4 : machineName_p="SetJyGridFT";
124 :
125 :
126 4 : }
127 :
128 24 : SetJyGridFT::~SetJyGridFT(){
129 :
130 24 : }
131 :
132 0 : SetJyGridFT::SetJyGridFT(const RecordInterface& stateRec)
133 0 : : GridFT()
134 : {
135 : // Construct from the input state record
136 0 : String error;
137 0 : if (!fromRecord(error, stateRec)) {
138 0 : throw (AipsError("Failed to create gridder: " + error));
139 : };
140 0 : }
141 8 : SetJyGridFT::SetJyGridFT(const SetJyGridFT& other) : GridFT()
142 : {
143 8 : machineName_p=("SetJyGridFT");
144 8 : operator=(other);
145 8 : }
146 :
147 8 : SetJyGridFT& SetJyGridFT::operator=(const SetJyGridFT& other)
148 : {
149 8 : if(this!=&other) {
150 : //Do the base parameters
151 8 : GridFT::operator=(other);
152 8 : freqscale_p.resize();
153 8 : freqscale_p=other.freqscale_p;
154 8 : scale_p.resize();
155 8 : scale_p=other.scale_p;
156 : }
157 8 : return *this;
158 : }
159 :
160 0 : FTMachine* SetJyGridFT::cloneFTM(){
161 0 : return new SetJyGridFT(*this);
162 : }
163 4 : void SetJyGridFT::setScale(const Vector<Double>& freq, const Vector<Double>& scale){
164 4 : freqscale_p.resize();
165 4 : freqscale_p=freq;
166 4 : scale_p.resize();
167 4 : scale_p=scale;
168 :
169 4 : }
170 644 : String SetJyGridFT::name() const{
171 :
172 644 : return String("SetJyGridFT");
173 : }
174 4 : void SetJyGridFT::initializeToVis(ImageInterface<Complex>& image,
175 : const VisBuffer& vb){
176 4 : setFreqInterpolation(String("nearest"));
177 4 : GridFT::initializeToVis(image, vb);
178 4 : }
179 :
180 0 : Bool SetJyGridFT::toRecord(String& error,
181 : RecordInterface& outRec, Bool withImage, const String diskimage)
182 : {
183 0 : if(!GridFT::toRecord(error, outRec, withImage, diskimage))
184 0 : return false;
185 0 : outRec.define("freqscale", freqscale_p);
186 0 : outRec.define("scaleamp", scale_p);
187 0 : return true;
188 : }
189 0 : Bool SetJyGridFT::fromRecord(String& error,
190 : const RecordInterface& inRec)
191 : {
192 :
193 0 : if(!GridFT::fromRecord(error, inRec))
194 0 : return false;
195 0 : freqscale_p.resize();
196 0 : inRec.get("freqscale", freqscale_p);
197 0 : scale_p.resize();
198 0 : inRec.get("scaleamp", scale_p);
199 0 : machineName_p="SetJyGridFT";
200 0 : return true;
201 : }
202 72 : void SetJyGridFT::get(VisBuffer& vb, Int row){
203 : //Did somebody really want this version.
204 72 : if(scale_p.nelements()==0)
205 0 : return GridFT::get(vb,row);
206 :
207 72 : gridOk(gridder->cSupport()(0));
208 : // If row is -1 then we pass through all rows
209 : Int startRow, endRow, nRow;
210 72 : if (row < 0) {
211 72 : nRow=vb.nRow();
212 72 : startRow=0;
213 72 : endRow=nRow-1;
214 : //unnecessary zeroing
215 : //vb.modelVisCube()=Complex(0.0,0.0);
216 : } else {
217 0 : nRow=1;
218 0 : startRow=row;
219 0 : endRow=row;
220 : //unnecessary zeroing
221 : //vb.modelVisCube().xyPlane(row)=Complex(0.0,0.0);
222 : }
223 :
224 : // Get the uvws in a form that Fortran can use
225 72 : Matrix<Double> uvw(3, vb.uvw().nelements());
226 72 : uvw=0.0;
227 72 : Vector<Double> dphase(vb.uvw().nelements());
228 72 : dphase=0.0;
229 : //NEGATING to correct for an image inversion problem
230 25344 : for (Int i=startRow;i<=endRow;i++) {
231 75816 : for (Int idim=0;idim<2;idim++) uvw(idim,i)=-vb.uvw()(i)(idim);
232 25272 : uvw(2,i)=vb.uvw()(i)(2);
233 : }
234 72 : rotateUVW(uvw, dphase, vb);
235 72 : refocus(uvw, vb.antenna1(), vb.antenna2(), dphase, vb);
236 :
237 :
238 :
239 : //Check if ms has changed then cache new spw and chan selection
240 72 : if(vb.newMS())
241 0 : matchAllSpwChans(vb);
242 :
243 :
244 : //Here we redo the match or use previous match
245 :
246 : //Channel matching for the actual spectral window of buffer
247 72 : if(doConversion_p[vb.spectralWindow()]){
248 0 : matchChannel(vb.spectralWindow(), vb);
249 : }
250 : else{
251 72 : chanMap.resize();
252 72 : chanMap=multiChanMap_p[vb.spectralWindow()];
253 : }
254 :
255 : //cerr << "chanMap " << chanMap << endl;
256 : //No point in reading data if its not matching in frequency
257 72 : if(max(chanMap)==-1)
258 0 : return;
259 :
260 144 : Cube<Complex> data;
261 144 : Cube<Int> flags;
262 72 : getInterpolateArrays(vb, data, flags);
263 :
264 144 : Vector<Double> lsrfreq;
265 : Bool convert;
266 72 : InterpolateArray1D<Double, Double>::InterpolationMethod meth= InterpolateArray1D<Double, Double>::nearestNeighbour;
267 72 : if(freqscale_p.nelements() > 2)
268 72 : meth= InterpolateArray1D<Double, Double>::linear;
269 72 : vb.lsrFrequency(vb.spectralWindow(), lsrfreq, convert);
270 72 : interpscale_p.resize(lsrfreq.nelements());
271 : InterpolateArray1D<Double,Double>::
272 72 : interpolate(interpscale_p,lsrfreq, freqscale_p, scale_p,meth);
273 :
274 :
275 :
276 : // IPosition s(data.shape());
277 72 : Int nvp=data.shape()(0);
278 72 : Int nvc=data.shape()(1);
279 72 : Int nvisrow=data.shape()(2);
280 :
281 :
282 : //cerr << "get flags " << min(flags) << " " << max(flags) << endl;
283 : Complex *datStorage;
284 : Bool isCopy;
285 72 : datStorage=data.getStorage(isCopy);
286 :
287 : ///
288 144 : Cube<Int> loc(2, nvc, nvisrow);
289 144 : Cube<Int> off(2, nvc, nRow);
290 144 : Matrix<Complex> phasor(nvc, nRow);
291 72 : Int csamp=gridder->cSampling();
292 : Bool delphase;
293 : Bool del;
294 72 : Complex * phasorstor=phasor.getStorage(delphase);
295 72 : const Double * visfreqstor=interpVisFreq_p.getStorage(del);
296 72 : const Double * scalestor=uvScale.getStorage(del);
297 72 : const Double * offsetstor=uvOffset.getStorage(del);
298 72 : const Double* uvstor= uvw.getStorage(del);
299 72 : Int * locstor=loc.getStorage(del);
300 72 : Int * offstor=off.getStorage(del);
301 72 : const Double *dpstor=dphase.getStorage(del);
302 72 : const Double *freqscalestor=interpscale_p.getStorage(del);
303 : //Vector<Double> f1=interpVisFreq_p.copy();
304 72 : Int nvchan=nvc;
305 : Int irow;
306 72 : #pragma omp parallel default(none) private(irow) firstprivate(visfreqstor, nvchan, scalestor, offsetstor, csamp, phasorstor, uvstor, locstor, offstor, dpstor) shared(startRow, endRow) num_threads(4)
307 : {
308 : #pragma omp for
309 : for (irow=startRow; irow<=endRow; ++irow){
310 : locateuvw(uvstor,dpstor, visfreqstor, nvchan, scalestor, offsetstor, csamp,
311 : locstor,
312 : offstor, phasorstor, irow);
313 : }
314 :
315 : }//end pragma parallel
316 :
317 72 : Int rbeg=startRow+1;
318 72 : Int rend=endRow+1;
319 :
320 :
321 144 : Vector<Int> rowFlags(vb.nRow());
322 72 : rowFlags=0;
323 72 : rowFlags(vb.flagRow())=true;
324 72 : if(!usezero_p) {
325 25344 : for (Int rownr=startRow; rownr<=endRow; rownr++) {
326 25272 : if(vb.antenna1()(rownr)==vb.antenna2()(rownr)) rowFlags(rownr)=1;
327 : }
328 : }
329 :
330 :
331 : //cerr <<"offset " << min(off) << " " <<max(off) << " length " << gridder->cFunction().shape() << endl;
332 :
333 : {
334 : Bool delgrid;
335 72 : const Complex* gridstor=griddedData.getStorage(delgrid);
336 72 : const Double * convfuncstor=(gridder->cFunction()).getStorage(del);
337 :
338 72 : const Int* pmapstor=polMap.getStorage(del);
339 72 : const Int *cmapstor=chanMap.getStorage(del);
340 72 : Int nc=nchan;
341 72 : Int np=npol;
342 72 : Int nxp=nx;
343 72 : Int nyp=ny;
344 72 : Int csupp=gridder->cSupport()(0);
345 72 : const Int * flagstor=flags.getStorage(del);
346 72 : const Int * rowflagstor=rowFlags.getStorage(del);
347 :
348 :
349 72 : Int npart=1;
350 : #ifdef _OPENMP
351 72 : Int nthreads=omp_get_max_threads();
352 72 : if (nthreads >3){
353 0 : npart=4;
354 : }
355 72 : else if(nthreads >1){
356 0 : npart=2;
357 : }
358 : #endif
359 72 : Int ix=0;
360 72 : #pragma omp parallel default(none) private(ix, rbeg, rend) firstprivate(datStorage, flagstor, rowflagstor, convfuncstor, pmapstor, cmapstor, gridstor, nxp, nyp, np, nc, csamp, csupp, nvp, nvc, nvisrow, phasorstor, locstor, offstor, freqscalestor) shared(npart) num_threads(npart)
361 : {
362 : #pragma omp for
363 : for (ix=0; ix< npart; ++ix){
364 : rbeg=ix*(nvisrow/npart)+1;
365 : rend=(ix != (npart-1)) ? (rbeg+(nvisrow/npart)-1) : (rbeg+(nvisrow/npart)+nvisrow%npart-1) ;
366 : //cerr << "rbeg " << rbeg << " rend " << rend << " " << nvisrow << endl;
367 :
368 : sectdgridjy(datStorage,
369 : &nvp,
370 : &nvc,
371 : flagstor,
372 : rowflagstor,
373 : &nvisrow,
374 : gridstor,
375 : &nxp,
376 : &nyp,
377 : &np,
378 : &nc,
379 : &csupp,
380 : &csamp,
381 : convfuncstor,
382 : cmapstor,
383 : pmapstor,
384 : &rbeg, &rend, locstor, offstor, phasorstor, freqscalestor);
385 : }//end pragma parallel
386 : }
387 72 : data.putStorage(datStorage, isCopy);
388 72 : griddedData.freeStorage(gridstor, delgrid);
389 : //cerr << "Get min max " << min(data) << " " << max(data) << endl;
390 : }
391 72 : interpolateFrequencyFromgrid(vb, data, FTMachine::MODEL);
392 :
393 :
394 : }
395 :
396 :
397 :
398 : }
|
