Line data Source code
1 : //# GridBoth.cc: Implementation of GridBoth 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/MeasurementComponents/GridBoth.h>
29 : #include <synthesis/TransformMachines/SimpCompGridMachine.h>
30 : #include <msvis/MSVis/VisBuffer.h>
31 : #include <casacore/casa/Arrays/ArrayMath.h>
32 : #include <casacore/images/Images/ImageInterface.h>
33 : #include <casacore/images/Images/TempImage.h>
34 :
35 : #include <synthesis/TransformMachines/StokesImageUtil.h>
36 :
37 : #include <components/ComponentModels/Flux.h>
38 : #include <components/ComponentModels/PointShape.h>
39 : #include <components/ComponentModels/ConstantSpectrum.h>
40 :
41 : #include <casacore/lattices/LRegions/LCBox.h>
42 : #include <casacore/lattices/LEL/LatticeExpr.h>
43 : #include <casacore/lattices/Lattices/SubLattice.h>
44 : #include <casacore/lattices/Lattices/LatticeIterator.h>
45 : #include <casacore/lattices/Lattices/LatticeStepper.h>
46 : #include <casacore/lattices/LEL/LatticeExpr.h>
47 : #include <casacore/casa/Containers/Record.h>
48 : #include <casacore/casa/BasicSL/String.h>
49 : #include <casacore/casa/Utilities/Assert.h>
50 : #include <casacore/casa/Exceptions/Error.h>
51 : #include <casacore/casa/OS/Timer.h>
52 : #include <sstream>
53 :
54 : using namespace casacore;
55 : namespace casa { //# NAMESPACE CASA - BEGIN
56 :
57 0 : GridBoth::GridBoth(SkyJones& sj, Long icachesize,
58 : Int itilesize,
59 : String sdConvType,
60 : String synConvType,
61 : Float padding,
62 : Float sdScale,
63 0 : Float sdWeight)
64 : : FTMachine(), synMachine_p(0), sdMachine_p(0), lastMachine_p(0),
65 0 : sdImage_p(0), synImage_p(0), sdScale_p(sdScale), sdWeight_p(sdWeight)
66 : {
67 0 : synMachine_p = new GridFT(icachesize, itilesize, synConvType,
68 0 : padding, false);
69 0 : sdMachine_p = new SDGrid(sj, icachesize, itilesize, sdConvType, -1);
70 0 : ok();
71 0 : }
72 :
73 0 : GridBoth::GridBoth(SkyJones& sj, Long icachesize,
74 : Int itilesize,
75 : MPosition mLocation,
76 : String sdConvType,
77 : String synConvType,
78 : Float padding,
79 : Float sdScale,
80 0 : Float sdWeight)
81 : : FTMachine(), synMachine_p(0), sdMachine_p(0), lastMachine_p(0),
82 0 : sdImage_p(0), synImage_p(0), sdScale_p(sdScale), sdWeight_p(sdWeight)
83 : {
84 0 : synMachine_p = new GridFT(icachesize, itilesize, synConvType, mLocation,
85 0 : padding, false);
86 0 : sdMachine_p = new SDGrid(mLocation, sj, icachesize, itilesize, sdConvType, -1);
87 0 : ok();
88 0 : }
89 :
90 0 : GridBoth::GridBoth(SkyJones& sj, Long icachesize,
91 : Int itilesize,
92 : MPosition mLocation,
93 : MDirection mDirection,
94 : String sdConvType,
95 : String synConvType,
96 : Float padding,
97 : Float sdScale,
98 0 : Float sdWeight)
99 :
100 : : FTMachine(), synMachine_p(0), sdMachine_p(0), lastMachine_p(0),
101 0 : sdImage_p(0), synImage_p(0), sdScale_p(sdScale), sdWeight_p(sdWeight)
102 : {
103 0 : synMachine_p = new GridFT(icachesize, itilesize, synConvType, mLocation,
104 0 : mDirection, padding, false);
105 0 : sdMachine_p = new SDGrid(mLocation, sj, icachesize, itilesize,
106 0 : sdConvType, -1);
107 0 : ok();
108 0 : }
109 :
110 0 : GridBoth::GridBoth(const RecordInterface& stateRec)
111 0 : : FTMachine()
112 : {
113 : // Construct from the input state record
114 0 : String error;
115 0 : if (!fromRecord(error, stateRec)) {
116 0 : throw (AipsError("Failed to create gridder: " + error));
117 : };
118 0 : }
119 :
120 : //----------------------------------------------------------------------
121 0 : GridBoth& GridBoth::operator=(const GridBoth& other)
122 : {
123 0 : if(this!=&other) {
124 0 : synMachine_p=other.synMachine_p;
125 0 : sdMachine_p=other.sdMachine_p;
126 0 : lastMachine_p=other.lastMachine_p;
127 0 : sdImage_p=other.sdImage_p;
128 0 : synImage_p=other.synImage_p;
129 0 : sdScale_p=other.sdScale_p;
130 0 : sdWeight_p=other.sdWeight_p;
131 : };
132 0 : return *this;
133 : };
134 :
135 : //----------------------------------------------------------------------
136 0 : GridBoth::GridBoth(const GridBoth& other):FTMachine()
137 : {
138 0 : operator=(other);
139 0 : }
140 :
141 0 : GridBoth::~GridBoth() {
142 0 : if(synMachine_p) delete synMachine_p; synMachine_p=0;
143 0 : if(sdMachine_p) delete sdMachine_p; sdMachine_p=0;
144 0 : if(sdImage_p) delete sdImage_p; sdImage_p=0;
145 0 : if(synImage_p) delete synImage_p; synImage_p=0;
146 0 : }
147 :
148 : //----------------------------------------------------------------------
149 0 : Bool GridBoth::changed(const VisBuffer& vb) {
150 0 : return synMachine_p->changed(vb)||sdMachine_p->changed(vb);
151 : }
152 :
153 : // Initialize for a transform from the Sky domain.
154 0 : void GridBoth::initializeToVis(ImageInterface<Complex>& iimage,
155 : const VisBuffer& vb)
156 : {
157 0 : ok();
158 :
159 0 : lastMachine_p=0;
160 :
161 0 : image=&iimage;
162 :
163 0 : if(sdImage_p) delete sdImage_p;
164 0 : sdImage_p=new TempImage<Complex>(iimage.shape(), iimage.coordinates());
165 0 : AlwaysAssert(sdImage_p, AipsError);
166 0 : sdImage_p->copyData(iimage);
167 0 : sdMachine_p->initializeToVis(*sdImage_p, vb);
168 :
169 0 : if(synImage_p) delete synImage_p;
170 0 : synImage_p=new TempImage<Complex>(iimage.shape(), iimage.coordinates());
171 0 : AlwaysAssert(synImage_p, AipsError);
172 0 : synImage_p->copyData(iimage);
173 0 : synMachine_p->initializeToVis(*synImage_p, vb);
174 :
175 0 : AlwaysAssert(sdImage_p->shape()==synImage_p->shape(), AipsError);
176 :
177 0 : }
178 :
179 0 : void GridBoth::finalizeToVis()
180 : {
181 0 : ok();
182 0 : synMachine_p->finalizeToVis();
183 0 : sdMachine_p->finalizeToVis();
184 0 : }
185 :
186 :
187 : // Initialize the transform to the Sky. Here we have to setup and initialize the
188 : // grid.
189 0 : void GridBoth::initializeToSky(ImageInterface<Complex>& iimage,
190 : Matrix<Float>& weight, const VisBuffer& vb)
191 : {
192 :
193 0 : ok();
194 :
195 0 : lastMachine_p=0;
196 :
197 0 : image=&iimage;
198 :
199 0 : Matrix<Float> sdWeights, synWeights;
200 :
201 0 : if(sdImage_p) delete sdImage_p;
202 :
203 :
204 0 : sdImage_p=new TempImage<Complex>(iimage.shape(), iimage.coordinates());
205 0 : AlwaysAssert(sdImage_p, AipsError);
206 0 : sdMachine_p->initializeToSky(*sdImage_p, sdWeights, vb);
207 :
208 0 : if(synImage_p) delete synImage_p;
209 0 : synImage_p=new TempImage<Complex>(iimage.shape(), iimage.coordinates());
210 0 : AlwaysAssert(synImage_p, AipsError);
211 0 : synMachine_p->initializeToSky(*synImage_p, synWeights, vb);
212 :
213 0 : sumWeight=0.0;
214 0 : weight.resize(synWeights.shape());
215 0 : weight=0.0;
216 :
217 :
218 :
219 0 : AlwaysAssert(sdImage_p->shape()==synImage_p->shape(), AipsError);
220 0 : }
221 :
222 0 : void GridBoth::finalizeToSky()
223 : {
224 0 : ok();
225 0 : synMachine_p->finalizeToSky();
226 0 : sdMachine_p->finalizeToSky();
227 0 : }
228 :
229 0 : void GridBoth::put(const VisBuffer& vb, Int row, Bool dopsf,
230 : FTMachine::Type type)
231 : {
232 0 : synMachine_p->put(vb, row, dopsf, type);
233 0 : sdMachine_p->put(vb, row, dopsf, type);
234 0 : }
235 :
236 0 : void GridBoth::get(VisBuffer& vb, Int row)
237 : {
238 0 : synMachine_p->get(vb, row);
239 0 : Cube<Complex> synModelVis(vb.modelVisCube().copy());
240 0 : sdMachine_p->get(vb, row);
241 0 : Cube<Complex> sdModelVis(vb.modelVisCube().copy());
242 0 : vb.setModelVisCube(sdModelVis+synModelVis);
243 0 : Cube<Complex> totalModelVis(vb.modelVisCube().copy());
244 0 : }
245 :
246 : // Finalize the transform to the Sky. We take the final images and
247 : // add them together, returning the resulting image
248 0 : ImageInterface<Complex>& GridBoth::getImage(Matrix<Float>& weights,
249 : Bool normalize)
250 : {
251 0 : ok();
252 0 : AlwaysAssert(image, AipsError);
253 0 : AlwaysAssert(synImage_p, AipsError);
254 0 : AlwaysAssert(sdImage_p, AipsError);
255 0 : AlwaysAssert(sdImage_p->shape()==synImage_p->shape(), AipsError);
256 :
257 0 : logIO() << LogOrigin("GridBoth", "getImage") << LogIO::NORMAL;
258 :
259 0 : Matrix<Float> synWeights, sdWeights;
260 :
261 0 : sdImage_p->copyData(sdMachine_p->getImage(sdWeights, false));
262 0 : synImage_p->copyData(synMachine_p->getImage(synWeights, false));
263 0 : Complex scale(sdScale_p*sdWeight_p);
264 0 : LatticeExpr<Complex> le(*synImage_p+scale*(*sdImage_p));
265 0 : image->copyData(le);
266 :
267 0 : if(normalize) {
268 0 : TempImage<Float> weightImage(image->shape(), image->coordinates());
269 0 : getWeightImage(weightImage, weights);
270 0 : if(max(weights)==0.0) {
271 : logIO() << LogIO::WARN << "No useful data in GridBoth: weights all zero"
272 0 : << LogIO::POST;
273 : }
274 0 : image->copyData((LatticeExpr<Complex>)(iif((weightImage<=0.0), 0.0, (*image)/(weightImage))));
275 : }
276 : else {
277 0 : weights.resize(synWeights.shape());
278 0 : weights=synWeights+sdWeight_p*sdWeights;
279 : }
280 :
281 0 : return *image;
282 : }
283 :
284 0 : void GridBoth::getWeightImage(ImageInterface<Float>& weightImage, Matrix<Float>& weights)
285 : {
286 0 : ok();
287 :
288 0 : logIO() << LogOrigin("GridBoth", "getWeightImage") << LogIO::NORMAL;
289 :
290 0 : AlwaysAssert(image, AipsError);
291 0 : AlwaysAssert(synImage_p, AipsError);
292 0 : AlwaysAssert(sdImage_p, AipsError);
293 0 : AlwaysAssert(sdImage_p->shape()==synImage_p->shape(), AipsError);
294 :
295 0 : logIO() << LogOrigin("GridBoth", "getWeightImage") << LogIO::NORMAL;
296 :
297 0 : TempImage<Float> sdWeightImage(sdImage_p->shape(), sdImage_p->coordinates());
298 0 : TempImage<Float> synWeightImage(synImage_p->shape(), synImage_p->coordinates());
299 :
300 0 : Matrix<Float> synWeights, sdWeights;
301 :
302 0 : sdMachine_p->getWeightImage(sdWeightImage, sdWeights);
303 0 : synMachine_p->getWeightImage(synWeightImage, synWeights);
304 :
305 0 : LatticeExpr<Float> le(synWeightImage+sdWeight_p*sdWeightImage);
306 0 : weightImage.copyData(le);
307 :
308 0 : weights.resize(synWeights.shape());
309 0 : weights=synWeights+sdWeight_p*sdWeights;
310 0 : }
311 :
312 : // Both these need filling out!
313 0 : Bool GridBoth::toRecord(String& error, RecordInterface& outRec,
314 : Bool withImage, const String diskimage ) {
315 0 : ok();
316 0 : String im1="";
317 0 : String im2="";
318 0 : if(diskimage != ""){
319 0 : im1=diskimage+"_synImage";
320 0 : im2=diskimage+"_sdImage";
321 : }
322 0 : return synMachine_p->toRecord(error, outRec, withImage, im1) &&
323 0 : sdMachine_p->toRecord(error, outRec, withImage, im2);
324 : }
325 :
326 0 : Bool GridBoth::fromRecord(String& error, const RecordInterface& inRec)
327 : {
328 0 : ok();
329 0 : return synMachine_p->fromRecord(error, inRec) &&
330 0 : sdMachine_p->fromRecord(error, inRec);
331 : }
332 :
333 0 : void GridBoth::ok() {
334 0 : AlwaysAssert(synMachine_p, AipsError);
335 0 : AlwaysAssert(sdMachine_p, AipsError);
336 0 : }
337 :
338 :
339 : } //# NAMESPACE CASA - END
340 :
|
