Line data Source code
1 : //# Mosaicft.cc: Implementation of MosaicFTNew class
2 : //# Copyright (C) 1997,1998,1999,2000,2001,2002,2003
3 : //# Associated Universities, Inc. Washington DC, USA.
4 : //#
5 : //# This library is free software; you can redistribute it and/or modify it
6 : //# under the terms of the GNU Library General Public License as published by
7 : //# the Free Software Foundation; either version 2 of the License, or (at your
8 : //# option) any later version.
9 : //#
10 : //# This library is distributed in the hope that it will be useful, but WITHOUT
11 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
13 : //# License for more details.
14 : //#
15 : //# You should have received a copy of the GNU Library General Public License
16 : //# along with this library; if not, write to the Free Software Foundation,
17 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
18 : //#
19 : //# Correspondence concerning AIPS++ should be addressed as follows:
20 : //# Internet email: casa-feedback@nrao.edu.
21 : //# Postal address: AIPS++ Project Office
22 : //# National Radio Astronomy Observatory
23 : //# 520 Edgemont Road
24 : //# Charlottesville, VA 22903-2475 USA
25 : //#
26 : //# $Id$
27 :
28 :
29 : #include <synthesis/TransformMachines2/MosaicFTNew.h>
30 : #include <msvis/MSVis/VisibilityIterator.h>
31 : #include <casacore/casa/Quanta/UnitMap.h>
32 : #include <casacore/casa/Quanta/MVTime.h>
33 : #include <casacore/casa/Quanta/UnitVal.h>
34 : #include <casacore/measures/Measures/Stokes.h>
35 : #include <casacore/measures/Measures/UVWMachine.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 :
45 : #include <synthesis/TransformMachines/StokesImageUtil.h>
46 :
47 : #include <casacore/images/Images/ImageInterface.h>
48 : #include <casacore/images/Images/PagedImage.h>
49 : #include <casacore/images/Images/SubImage.h>
50 : #include <casacore/images/Regions/ImageRegion.h>
51 : #include <casacore/images/Regions/WCBox.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/Slice.h>
60 : #include <casacore/casa/Arrays/Matrix.h>
61 : #include <casacore/casa/Arrays/Cube.h>
62 : #include <casacore/casa/Arrays/MatrixIter.h>
63 : #include <casacore/casa/BasicSL/String.h>
64 : #include <casacore/casa/Utilities/Assert.h>
65 : #include <casacore/casa/Exceptions/Error.h>
66 : #include <casacore/lattices/Lattices/ArrayLattice.h>
67 : #include <casacore/lattices/Lattices/SubLattice.h>
68 : #include <casacore/lattices/LRegions/LCBox.h>
69 : #include <casacore/lattices/LEL/LatticeExpr.h>
70 : #include <casacore/lattices/Lattices/LatticeCache.h>
71 : #include <casacore/lattices/LatticeMath/LatticeFFT.h>
72 : #include <casacore/lattices/Lattices/LatticeIterator.h>
73 : #include <casacore/lattices/Lattices/LatticeStepper.h>
74 : #include <casacore/casa/Utilities/CompositeNumber.h>
75 : #include <casacore/casa/OS/Timer.h>
76 : #include <casacore/casa/OS/HostInfo.h>
77 : #include <sstream>
78 : #ifdef _OPENMP
79 : #include <omp.h>
80 : #endif
81 : using namespace casacore;
82 : namespace casa { //# NAMESPACE CASA - BEGIN
83 : namespace refim {//# namespace for imaging refactor
84 : using namespace casacore;
85 : using namespace casa;
86 : using namespace casacore;
87 : using namespace casa::refim;
88 :
89 69 : refim::FTMachine* MosaicFTNew::cloneFTM(){
90 69 : return new MosaicFTNew(*this);
91 : }
92 :
93 0 : MosaicFTNew::MosaicFTNew(const RecordInterface& stateRec)
94 0 : : MosaicFT(stateRec)
95 : {
96 :
97 0 : }
98 0 : MosaicFTNew& MosaicFTNew::operator=(const MosaicFTNew& other) {
99 0 : if(this!=&other) {
100 :
101 : //Do the base parameters
102 0 : MosaicFT::operator=(other);
103 :
104 : }
105 0 : return *this;
106 :
107 : }
108 : // Finalize the FFT to the Sky. Here we actually do the FFT and
109 : // return the resulting image
110 289 : ImageInterface<Complex>& MosaicFTNew::getImage(Matrix<Float>& weights,
111 : Bool normalize)
112 : {
113 : //AlwaysAssert(lattice, AipsError);
114 289 : AlwaysAssert(image, AipsError);
115 :
116 289 : if((griddedData.nelements() ==0) && (griddedData2.nelements()==0)){
117 0 : throw(AipsError("Programmer error ...request for image without right order of calls"));
118 : }
119 :
120 289 : logIO() << LogOrigin("MosaicFTNew", "getImage") << LogIO::NORMAL;
121 :
122 289 : weights.resize(sumWeight.shape());
123 :
124 289 : convertArray(weights, sumWeight);
125 :
126 : // If the weights are all zero then we cannot normalize
127 : // otherwise we don't care.
128 289 : if(max(weights)==0.0) {
129 0 : if(normalize) {
130 : logIO() << LogIO::SEVERE << "No useful data in MosaicFTNew: weights all zero"
131 0 : << LogIO::POST;
132 : }
133 : else {
134 : logIO() << LogIO::WARN << "No useful data in MosaicFTNew: weights all zero"
135 0 : << LogIO::POST;
136 0 : image->set(0.0);
137 : }
138 : }
139 : else {
140 :
141 : //const IPosition latticeShape = lattice->shape();
142 :
143 : logIO() << LogIO::DEBUGGING
144 289 : << "Starting FFT and scaling of image" << LogIO::POST;
145 289 : if(useDoubleGrid_p){
146 289 : ArrayLattice<DComplex> darrayLattice(griddedData2);
147 289 : ft_p.c2cFFT(darrayLattice,false);
148 289 : griddedData.resize(griddedData2.shape());
149 289 : convertArray(griddedData, griddedData2);
150 :
151 : //Don't need the double-prec grid anymore...
152 289 : griddedData2.resize();
153 289 : arrayLattice = new ArrayLattice<Complex>(griddedData);
154 289 : lattice=arrayLattice;
155 :
156 289 : }
157 : else{
158 0 : arrayLattice = new ArrayLattice<Complex>(griddedData);
159 0 : lattice=arrayLattice;
160 0 : ft_p.c2cFFT(*lattice,false);
161 : }
162 :
163 : {
164 289 : Int inx = lattice->shape()(0);
165 289 : Int iny = lattice->shape()(1);
166 289 : Vector<Complex> correction(inx);
167 289 : correction=Complex(1.0, 0.0);
168 289 : Vector<Float> sincConvX(nx);
169 262469 : for (Int ix=0;ix<nx;ix++) {
170 262180 : Float x=C::pi*Float(ix-nx/2)/(Float(nx)*Float(convSampling));
171 262180 : if(ix==nx/2) {
172 289 : sincConvX(ix)=1.0;
173 : }
174 : else {
175 261891 : sincConvX(ix)=sin(x)/x;
176 : }
177 : }
178 289 : Vector<Float> sincConvY(ny);
179 261989 : for (Int ix=0;ix<ny;ix++) {
180 261700 : Float x=C::pi*Float(ix-ny/2)/(Float(ny)*Float(convSampling));
181 261700 : if(ix==ny/2) {
182 289 : sincConvY(ix)=1.0;
183 : }
184 : else {
185 261411 : sincConvY(ix)=sin(x)/x;
186 : }
187 : }
188 :
189 :
190 289 : IPosition cursorShape(4, inx, 1, 1, 1);
191 289 : IPosition axisPath(4, 0, 1, 2, 3);
192 289 : LatticeStepper lsx(lattice->shape(), cursorShape, axisPath);
193 289 : LatticeIterator<Complex> lix(*lattice, lsx);
194 497565 : for(lix.reset();!lix.atEnd();lix++)
195 : {
196 497276 : Int pol=lix.position()(2);
197 497276 : Int chan=lix.position()(3);
198 497276 : Int iy=lix.position()(1);
199 541497772 : for (Int ix=0;ix<nx;ix++) {
200 541000496 : correction(ix)=1.0/(sincConvX(ix)*sincConvY(iy));
201 : }
202 497276 : if(normalize)
203 : {
204 2048 : if(weights(pol, chan)!=0.0)
205 : {
206 2048 : Complex rnorm(Float(inx)*Float(iny)/weights(pol,chan));
207 2048 : lix.rwVectorCursor()*=rnorm*correction;
208 : }
209 : else {
210 0 : lix.woCursor()=0.0;
211 : }
212 : } else {
213 495228 : Complex rnorm(Float(inx)*Float(iny));
214 495228 : lix.rwVectorCursor()*=rnorm*correction;
215 : }
216 : }
217 289 : }
218 :
219 : //if(!isTiled)
220 : {
221 289 : LatticeLocker lock1 (*(image), FileLocker::Write);
222 : // Check the section from the image BEFORE converting to a lattice
223 578 : IPosition blc(4, (nx-image->shape()(0)+(nx%2==0))/2,
224 578 : (ny-image->shape()(1)+(ny%2==0))/2, 0, 0);
225 289 : IPosition stride(4, 1);
226 578 : IPosition trc(blc+image->shape()-stride);
227 :
228 : // Do the copy
229 289 : IPosition start(4, 0);
230 289 : image->put(griddedData(blc, trc));
231 289 : }
232 : }
233 289 : if(!arrayLattice.null()) arrayLattice=0;
234 289 : if(!lattice.null()) lattice=0;
235 289 : griddedData.resize();
236 289 : image->clearCache();
237 289 : return *image;
238 : }
239 :
240 : // Get weight image
241 113 : void MosaicFTNew::getWeightImage(ImageInterface<Float>& weightImage,
242 : Matrix<Float>& weights)
243 : {
244 :
245 113 : logIO() << LogOrigin("MosaicFTNew", "getWeightImage") << LogIO::NORMAL;
246 :
247 : //cerr << "SUMWEIGHT " << sumWeight << endl;
248 113 : weights.resize(sumWeight.shape());
249 113 : convertArray(weights,sumWeight);
250 113 : if(!skyCoverage_p)
251 0 : throw(AipsError("MosaicFTNew::getWeightImage: called before initializing"));
252 : //cerr << "skyCoverage_p " << skyCoverage_p.get() << endl;
253 113 : Record rec=skyCoverage_p->miscInfo();
254 113 : Float inx=1;
255 113 : Float iny=1;
256 113 : Bool isscaled=rec.isDefined("isscaled") && rec.asBool("isscaled");
257 : //cerr << "isscaled " << isscaled << " max " << max(skyCoverage_p->get()) << endl;
258 113 : if( !isscaled){
259 113 : inx = skyCoverage_p->shape()(0);
260 113 : iny = skyCoverage_p->shape()(1);
261 : }
262 : //cerr << "inx, iny " << inx << " " << iny << endl;
263 113 : weightImage.copyData( (LatticeExpr<Float>) ( (*skyCoverage_p)*inx*iny ) );
264 :
265 113 : skyCoverage_p->tempClose();
266 :
267 113 : }
268 :
269 :
270 : } // REFIM ends
271 : } //# NAMESPACE CASA - END
|
