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: 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 : //# $Id$
27 :
28 : #include <synthesis/TransformMachines2/MosaicFTNew.h>
29 :
30 :
31 : #include <msvis/MSVis/VisibilityIterator.h>
32 : #include <casacore/casa/Quanta/UnitMap.h>
33 : #include <casacore/casa/Quanta/MVTime.h>
34 : #include <casacore/casa/Quanta/UnitVal.h>
35 : #include <casacore/measures/Measures/Stokes.h>
36 : #include <casacore/measures/Measures/UVWMachine.h>
37 : #include <casacore/coordinates/Coordinates/CoordinateSystem.h>
38 : #include <casacore/coordinates/Coordinates/DirectionCoordinate.h>
39 : #include <casacore/coordinates/Coordinates/SpectralCoordinate.h>
40 : #include <casacore/coordinates/Coordinates/StokesCoordinate.h>
41 : #include <casacore/coordinates/Coordinates/Projection.h>
42 : #include <casacore/ms/MeasurementSets/MSColumns.h>
43 : #include <casacore/casa/BasicSL/Constants.h>
44 : #include <casacore/scimath/Mathematics/FFTServer.h>
45 :
46 : #include <synthesis/TransformMachines/StokesImageUtil.h>
47 :
48 : #include <casacore/images/Images/ImageInterface.h>
49 : #include <casacore/images/Images/PagedImage.h>
50 : #include <casacore/images/Images/SubImage.h>
51 : #include <casacore/images/Regions/ImageRegion.h>
52 : #include <casacore/images/Regions/WCBox.h>
53 : #include <casacore/casa/Containers/Block.h>
54 : #include <casacore/casa/Containers/Record.h>
55 : #include <casacore/casa/Arrays/ArrayLogical.h>
56 : #include <casacore/casa/Arrays/ArrayMath.h>
57 : #include <casacore/casa/Arrays/Array.h>
58 : #include <casacore/casa/Arrays/MaskedArray.h>
59 : #include <casacore/casa/Arrays/Vector.h>
60 : #include <casacore/casa/Arrays/Slice.h>
61 : #include <casacore/casa/Arrays/Matrix.h>
62 : #include <casacore/casa/Arrays/Cube.h>
63 : #include <casacore/casa/Arrays/MatrixIter.h>
64 : #include <casacore/casa/BasicSL/String.h>
65 : #include <casacore/casa/Utilities/Assert.h>
66 : #include <casacore/casa/Exceptions/Error.h>
67 : #include <casacore/lattices/Lattices/ArrayLattice.h>
68 : #include <casacore/lattices/Lattices/SubLattice.h>
69 : #include <casacore/lattices/LRegions/LCBox.h>
70 : #include <casacore/lattices/LEL/LatticeExpr.h>
71 : #include <casacore/lattices/Lattices/LatticeCache.h>
72 : #include <casacore/lattices/LatticeMath/LatticeFFT.h>
73 : #include <casacore/lattices/Lattices/LatticeIterator.h>
74 : #include <casacore/lattices/Lattices/LatticeStepper.h>
75 : #include <casacore/casa/Utilities/CompositeNumber.h>
76 : #include <casacore/casa/OS/Timer.h>
77 : #include <casacore/casa/OS/HostInfo.h>
78 : #include <sstream>
79 : #ifdef _OPENMP
80 : #include <omp.h>
81 : #endif
82 : using namespace casacore;
83 : namespace casa { //# NAMESPACE CASA - BEGIN
84 : namespace refim {//# namespace for imaging refactor
85 : using namespace casacore;
86 : using namespace casa;
87 : using namespace casacore;
88 : using namespace casa::refim;
89 :
90 69 : refim::FTMachine* MosaicFTNew::cloneFTM(){
91 69 : return new MosaicFTNew(*this);
92 : }
93 :
94 0 : MosaicFTNew::MosaicFTNew(const RecordInterface& stateRec)
95 0 : : MosaicFT(stateRec)
96 : {
97 :
98 0 : }
99 : // Finalize the FFT to the Sky. Here we actually do the FFT and
100 : // return the resulting image
101 276 : ImageInterface<Complex>& MosaicFTNew::getImage(Matrix<Float>& weights,
102 : Bool normalize)
103 : {
104 : //AlwaysAssert(lattice, AipsError);
105 276 : AlwaysAssert(image, AipsError);
106 :
107 276 : if((griddedData.nelements() ==0) && (griddedData2.nelements()==0)){
108 0 : throw(AipsError("Programmer error ...request for image without right order of calls"));
109 : }
110 :
111 276 : logIO() << LogOrigin("MosaicFTNew", "getImage") << LogIO::NORMAL;
112 :
113 276 : weights.resize(sumWeight.shape());
114 :
115 276 : convertArray(weights, sumWeight);
116 :
117 : // If the weights are all zero then we cannot normalize
118 : // otherwise we don't care.
119 276 : if(max(weights)==0.0) {
120 0 : if(normalize) {
121 : logIO() << LogIO::SEVERE << "No useful data in MosaicFTNew: weights all zero"
122 0 : << LogIO::POST;
123 : }
124 : else {
125 : logIO() << LogIO::WARN << "No useful data in MosaicFTNew: weights all zero"
126 0 : << LogIO::POST;
127 0 : image->set(0.0);
128 : }
129 : }
130 : else {
131 :
132 : //const IPosition latticeShape = lattice->shape();
133 :
134 : logIO() << LogIO::DEBUGGING
135 276 : << "Starting FFT and scaling of image" << LogIO::POST;
136 276 : if(useDoubleGrid_p){
137 552 : ArrayLattice<DComplex> darrayLattice(griddedData2);
138 276 : ft_p.c2cFFT(darrayLattice,false);
139 276 : griddedData.resize(griddedData2.shape());
140 276 : convertArray(griddedData, griddedData2);
141 :
142 : //Don't need the double-prec grid anymore...
143 276 : griddedData2.resize();
144 276 : arrayLattice = new ArrayLattice<Complex>(griddedData);
145 276 : lattice=arrayLattice;
146 :
147 : }
148 : else{
149 0 : arrayLattice = new ArrayLattice<Complex>(griddedData);
150 0 : lattice=arrayLattice;
151 0 : ft_p.c2cFFT(*lattice,false);
152 : }
153 :
154 : {
155 276 : Int inx = lattice->shape()(0);
156 276 : Int iny = lattice->shape()(1);
157 552 : Vector<Complex> correction(inx);
158 276 : correction=Complex(1.0, 0.0);
159 552 : Vector<Float> sincConvX(nx);
160 255928 : for (Int ix=0;ix<nx;ix++) {
161 255652 : Float x=C::pi*Float(ix-nx/2)/(Float(nx)*Float(convSampling));
162 255652 : if(ix==nx/2) {
163 276 : sincConvX(ix)=1.0;
164 : }
165 : else {
166 255376 : sincConvX(ix)=sin(x)/x;
167 : }
168 : }
169 552 : Vector<Float> sincConvY(ny);
170 255688 : for (Int ix=0;ix<ny;ix++) {
171 255412 : Float x=C::pi*Float(ix-ny/2)/(Float(ny)*Float(convSampling));
172 255412 : if(ix==ny/2) {
173 276 : sincConvY(ix)=1.0;
174 : }
175 : else {
176 255136 : sincConvY(ix)=sin(x)/x;
177 : }
178 : }
179 :
180 :
181 552 : IPosition cursorShape(4, inx, 1, 1, 1);
182 552 : IPosition axisPath(4, 0, 1, 2, 3);
183 552 : LatticeStepper lsx(lattice->shape(), cursorShape, axisPath);
184 552 : LatticeIterator<Complex> lix(*lattice, lsx);
185 480632 : for(lix.reset();!lix.atEnd();lix++)
186 : {
187 480356 : Int pol=lix.position()(2);
188 480356 : Int chan=lix.position()(3);
189 480356 : Int iy=lix.position()(1);
190 533113492 : for (Int ix=0;ix<nx;ix++) {
191 532633136 : correction(ix)=1.0/(sincConvX(ix)*sincConvY(iy));
192 : }
193 480356 : if(normalize)
194 : {
195 2048 : if(weights(pol, chan)!=0.0)
196 : {
197 2048 : Complex rnorm(Float(inx)*Float(iny)/weights(pol,chan));
198 2048 : lix.rwVectorCursor()*=rnorm*correction;
199 : }
200 : else {
201 0 : lix.woCursor()=0.0;
202 : }
203 : } else {
204 478308 : Complex rnorm(Float(inx)*Float(iny));
205 478308 : lix.rwVectorCursor()*=rnorm*correction;
206 : }
207 : }
208 : }
209 :
210 : //if(!isTiled)
211 : {
212 552 : LatticeLocker lock1 (*(image), FileLocker::Write);
213 : // Check the section from the image BEFORE converting to a lattice
214 552 : IPosition blc(4, (nx-image->shape()(0)+(nx%2==0))/2,
215 828 : (ny-image->shape()(1)+(ny%2==0))/2, 0, 0);
216 552 : IPosition stride(4, 1);
217 828 : IPosition trc(blc+image->shape()-stride);
218 :
219 : // Do the copy
220 276 : IPosition start(4, 0);
221 276 : image->put(griddedData(blc, trc));
222 : }
223 : }
224 276 : if(!arrayLattice.null()) arrayLattice=0;
225 276 : if(!lattice.null()) lattice=0;
226 276 : griddedData.resize();
227 276 : image->clearCache();
228 276 : return *image;
229 : }
230 :
231 : // Get weight image
232 107 : void MosaicFTNew::getWeightImage(ImageInterface<Float>& weightImage,
233 : Matrix<Float>& weights)
234 : {
235 :
236 107 : logIO() << LogOrigin("MosaicFTNew", "getWeightImage") << LogIO::NORMAL;
237 :
238 : //cerr << "SUMWEIGHT " << sumWeight << endl;
239 107 : weights.resize(sumWeight.shape());
240 107 : convertArray(weights,sumWeight);
241 107 : if(!skyCoverage_p)
242 0 : throw(AipsError("MosaicFTNew::getWeightImage: called before initializing"));
243 : //cerr << "skyCoverage_p " << skyCoverage_p.get() << endl;
244 214 : Record rec=skyCoverage_p->miscInfo();
245 107 : Float inx=1;
246 107 : Float iny=1;
247 107 : Bool isscaled=rec.isDefined("isscaled") && rec.asBool("isscaled");
248 : //cerr << "isscaled " << isscaled << " max " << max(skyCoverage_p->get()) << endl;
249 107 : if( !isscaled){
250 107 : inx = skyCoverage_p->shape()(0);
251 107 : iny = skyCoverage_p->shape()(1);
252 : }
253 : //cerr << "inx, iny " << inx << " " << iny << endl;
254 107 : weightImage.copyData( (LatticeExpr<Float>) ( (*skyCoverage_p)*inx*iny ) );
255 :
256 107 : skyCoverage_p->tempClose();
257 :
258 107 : }
259 :
260 :
261 : } // REFIM ends
262 : } //# NAMESPACE CASA - END
|
