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 287 : ImageInterface<Complex>& MosaicFTNew::getImage(Matrix<Float>& weights,
111 : Bool normalize)
112 : {
113 : //AlwaysAssert(lattice, AipsError);
114 287 : AlwaysAssert(image, AipsError);
115 :
116 287 : if((griddedData.nelements() ==0) && (griddedData2.nelements()==0)){
117 0 : throw(AipsError("Programmer error ...request for image without right order of calls"));
118 : }
119 :
120 287 : logIO() << LogOrigin("MosaicFTNew", "getImage") << LogIO::NORMAL;
121 :
122 287 : weights.resize(sumWeight.shape());
123 :
124 287 : convertArray(weights, sumWeight);
125 :
126 : // If the weights are all zero then we cannot normalize
127 : // otherwise we don't care.
128 287 : 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 287 : << "Starting FFT and scaling of image" << LogIO::POST;
145 287 : if(useDoubleGrid_p){
146 287 : ArrayLattice<DComplex> darrayLattice(griddedData2);
147 287 : ft_p.c2cFFT(darrayLattice,false);
148 287 : griddedData.resize(griddedData2.shape());
149 287 : convertArray(griddedData, griddedData2);
150 :
151 : //Don't need the double-prec grid anymore...
152 287 : griddedData2.resize();
153 287 : arrayLattice = new ArrayLattice<Complex>(griddedData);
154 287 : lattice=arrayLattice;
155 :
156 287 : }
157 : else{
158 0 : arrayLattice = new ArrayLattice<Complex>(griddedData);
159 0 : lattice=arrayLattice;
160 0 : ft_p.c2cFFT(*lattice,false);
161 : }
162 :
163 : {
164 287 : Int inx = lattice->shape()(0);
165 287 : Int iny = lattice->shape()(1);
166 287 : Vector<Complex> correction(inx);
167 287 : correction=Complex(1.0, 0.0);
168 287 : Vector<Float> sincConvX(nx);
169 261443 : for (Int ix=0;ix<nx;ix++) {
170 261156 : Float x=C::pi*Float(ix-nx/2)/(Float(nx)*Float(convSampling));
171 261156 : if(ix==nx/2) {
172 287 : sincConvX(ix)=1.0;
173 : }
174 : else {
175 260869 : sincConvX(ix)=sin(x)/x;
176 : }
177 : }
178 287 : Vector<Float> sincConvY(ny);
179 260963 : for (Int ix=0;ix<ny;ix++) {
180 260676 : Float x=C::pi*Float(ix-ny/2)/(Float(ny)*Float(convSampling));
181 260676 : if(ix==ny/2) {
182 287 : sincConvY(ix)=1.0;
183 : }
184 : else {
185 260389 : sincConvY(ix)=sin(x)/x;
186 : }
187 : }
188 :
189 :
190 287 : IPosition cursorShape(4, inx, 1, 1, 1);
191 287 : IPosition axisPath(4, 0, 1, 2, 3);
192 287 : LatticeStepper lsx(lattice->shape(), cursorShape, axisPath);
193 287 : LatticeIterator<Complex> lix(*lattice, lsx);
194 496539 : for(lix.reset();!lix.atEnd();lix++)
195 : {
196 496252 : Int pol=lix.position()(2);
197 496252 : Int chan=lix.position()(3);
198 496252 : Int iy=lix.position()(1);
199 540972460 : for (Int ix=0;ix<nx;ix++) {
200 540476208 : correction(ix)=1.0/(sincConvX(ix)*sincConvY(iy));
201 : }
202 496252 : 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 494204 : Complex rnorm(Float(inx)*Float(iny));
214 494204 : lix.rwVectorCursor()*=rnorm*correction;
215 : }
216 : }
217 287 : }
218 :
219 : //if(!isTiled)
220 : {
221 287 : LatticeLocker lock1 (*(image), FileLocker::Write);
222 : // Check the section from the image BEFORE converting to a lattice
223 574 : IPosition blc(4, (nx-image->shape()(0)+(nx%2==0))/2,
224 574 : (ny-image->shape()(1)+(ny%2==0))/2, 0, 0);
225 287 : IPosition stride(4, 1);
226 574 : IPosition trc(blc+image->shape()-stride);
227 :
228 : // Do the copy
229 287 : IPosition start(4, 0);
230 287 : image->put(griddedData(blc, trc));
231 287 : }
232 : }
233 287 : if(!arrayLattice.null()) arrayLattice=0;
234 287 : if(!lattice.null()) lattice=0;
235 287 : griddedData.resize();
236 287 : image->clearCache();
237 287 : return *image;
238 : }
239 :
240 : // Get weight image
241 112 : void MosaicFTNew::getWeightImage(ImageInterface<Float>& weightImage,
242 : Matrix<Float>& weights)
243 : {
244 :
245 112 : logIO() << LogOrigin("MosaicFTNew", "getWeightImage") << LogIO::NORMAL;
246 :
247 : //cerr << "SUMWEIGHT " << sumWeight << endl;
248 112 : weights.resize(sumWeight.shape());
249 112 : convertArray(weights,sumWeight);
250 112 : if(!skyCoverage_p)
251 0 : throw(AipsError("MosaicFTNew::getWeightImage: called before initializing"));
252 : //cerr << "skyCoverage_p " << skyCoverage_p.get() << endl;
253 112 : Record rec=skyCoverage_p->miscInfo();
254 112 : Float inx=1;
255 112 : Float iny=1;
256 112 : Bool isscaled=rec.isDefined("isscaled") && rec.asBool("isscaled");
257 : //cerr << "isscaled " << isscaled << " max " << max(skyCoverage_p->get()) << endl;
258 112 : if( !isscaled){
259 112 : inx = skyCoverage_p->shape()(0);
260 112 : iny = skyCoverage_p->shape()(1);
261 : }
262 : //cerr << "inx, iny " << inx << " " << iny << endl;
263 112 : weightImage.copyData( (LatticeExpr<Float>) ( (*skyCoverage_p)*inx*iny ) );
264 :
265 112 : skyCoverage_p->tempClose();
266 :
267 112 : }
268 :
269 :
270 : } // REFIM ends
271 : } //# NAMESPACE CASA - END
|
