Line data Source code
1 : //# SkyEquation.h: SkyEquation definition
2 : //# Copyright (C) 1996,1997,1998,1999,2000,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 adressed 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 : //#
27 : //# $Id$
28 :
29 : #ifndef SYNTHESIS_SKYEQUATION_H
30 : #define SYNTHESIS_SKYEQUATION_H
31 :
32 : #include <casacore/casa/aips.h>
33 : #include <casacore/casa/Arrays/Matrix.h>
34 : #include <casacore/casa/BasicSL/Complex.h>
35 : #include <casacore/casa/Logging/LogMessage.h>
36 : #include <casacore/casa/Logging/LogSink.h>
37 : #include <synthesis/TransformMachines/FTMachine.h>
38 : #include <msvis/MSVis/VisBuffer.h>
39 : #include <casacore/ms/MeasurementSets/MSMainEnums.h>
40 :
41 : namespace casacore{
42 :
43 : class UVWMachine;
44 : template <class T> class ImageInterface;
45 : template <class T> class TempImage;
46 : }
47 :
48 : namespace casa { //# NAMESPACE CASA - BEGIN
49 :
50 : // <summary> Relate Sky brightness to the visibility </summary>
51 :
52 : // <use visibility=export>
53 :
54 : // <reviewed reviewer="" date="" tests="" demos="">
55 :
56 : // <prerequisite>
57 : // <li> <linkto module="MeasurementComponents">MeasurementComponents</linkto> module
58 : // </prerequisite>
59 : //
60 : // <etymology>
61 : // Sky Equation encapsulates the equation between the sky brightness
62 : // and the visibility (or coherence) measured by a generic instrument
63 : // </etymology>
64 : //
65 : // <synopsis>
66 : // This is responsible for the Sky-based part of Measurement Equation of the Generic
67 : // Interferometer due to Hamaker, Bregman and Sault and later extended
68 : // by Noordam, and Cornwell.
69 : //
70 : // See <linkto module="MeasurementEquations">MeasurementEquations</linkto>
71 : // for more details of the form of the SkyEquation.
72 : //
73 : // The principal use of SkyEquation is that, as described in
74 : // <linkto module="MeasurementEquations">MeasurementEquations</linkto>,
75 : // the gradients of chi-squared may be calculated and returned
76 : // as an image.
77 : //
78 : // The following components can be plugged into SkyEquation
79 : // <ul>
80 : // <li> Antenna-based direction-dependent terms: <linkto class="SkyJones">SkyJones</linkto>
81 : // <li> Sky brightness model: <linkto class="SkyModel">SkyModel</linkto>
82 : // <li> Fourier transform machine: <linkto class="FTMachine">FTMachine</linkto>
83 : // </ul>
84 : // </synopsis>
85 : //
86 : // <example>
87 : // <srcblock>
88 : //
89 : // // Read the VisSet from disk
90 : // VisSet vs("3c84.MS");
91 : //
92 : // casacore::PagedImage<casacore::Float> im("3c84.modelImage");
93 : //
94 : // // Create an ImageSkyModel from an image on disk
95 : // ImageSkyModel ism(im);
96 : //
97 : // // FTMachine
98 : // GridFT ft;
99 : //
100 : // SkyEquation se(ism, vs, ft);
101 : //
102 : // // Make a Clean Image and write it out
103 : // HogbomCleanImageSkyModel csm(ism);
104 : // if (csm.solveSkyModel()) {
105 : // casacore::PagedImage<casacore::Float> cleanImage=csm.getImage();
106 : // cleanImage.setName("3c84.cleanImage");
107 : // }
108 : //
109 : // </srcblock>
110 : // </example>
111 : //
112 : // <motivation>
113 : // SkyEquation is needed to encapsulate part of the measurement equation for
114 : // both single dish and synthesis observations. The idea is that the structure of many
115 : // imaging algorithms is much the same for many different types of telescope.
116 : // SkyEquation is part of a framework of classes that are
117 : // designed for synthesis and single dish imaging. The others are the
118 : // <linkto module=MeasurementComponents>MeasurementComponents</linkto>.
119 : //
120 : // </motivation>
121 : //
122 : // <todo asof="98/2/17">
123 : // <li> Implement SkyJones
124 : // </todo>
125 :
126 : // Forward declarations
127 : class SkyJones;
128 : class SkyModel;
129 : class VisSet;
130 : class FTMachine;
131 : class SkyComponent;
132 : class ComponentList;
133 : class ComponentFTMachine;
134 : class ROVisibilityIterator;
135 : class VisibilityIterator;
136 :
137 :
138 : class SkyEquation {
139 : public:
140 :
141 : // Define a SkyEquation linking a VisSet vs with a SkyModel sm
142 : // using an FTMachine ft for transforms in both directions
143 : SkyEquation(SkyModel& sm, VisSet& vs, FTMachine& ft);
144 :
145 : // Define a SkyEquation linking a VisSet vs with a SkyModel sm
146 : // using an FTMachine ft for transforms in both directions and
147 : // a ComponentFTMachine for the component lists
148 : SkyEquation(SkyModel& sm, VisSet& vs, FTMachine& ft,
149 : ComponentFTMachine& cft, casacore::Bool noModelcol=false);
150 :
151 :
152 : //SkyEquation with ROVisIter
153 : SkyEquation(SkyModel& sm, ROVisibilityIterator& vi, FTMachine& ft,
154 : ComponentFTMachine& cft, casacore::Bool noModelCol);
155 :
156 : // Define a SkyEquation linking a VisSet vs with a SkyModel sm
157 : // using an FTMachine ft for Sky->Vis and ift for Vis->Sky
158 : SkyEquation(SkyModel& sm, VisSet& vs, FTMachine& ft, FTMachine& ift);
159 :
160 : // Define a SkyEquation linking a VisSet vs with a SkyModel sm
161 : // using an FTMachine ft for Sky->Vis and ift for Vis->Sky and
162 : // a ComponentFTMachine for the component lists
163 : // Default make use of MODEL_DATA Column of ms
164 : SkyEquation(SkyModel& sm, VisSet& vs, FTMachine& ft, FTMachine& ift,
165 : ComponentFTMachine& cft);
166 :
167 :
168 : // Return the SkyModel used by this SkyEquation
169 : SkyModel& skyModel() {return *sm_;};
170 :
171 : // Return the VisSet used by this SkyEquation
172 : VisSet& visSet() {return *vs_;};
173 :
174 : // Return the (Sky->Vis) FTMachine used by this SkyEquation
175 0 : FTMachine& fTMachine() {return *ft_;};
176 :
177 : // Return the (Vis->Sky) FTMachine used by this SkyEquation
178 : FTMachine& iFTMachine() {return *ift_;};
179 :
180 : // Return the (Sky->Vis) ComponentFTMachine used by this SkyEquation
181 : ComponentFTMachine& componentFTMachine() {return *cft_;};
182 :
183 : // Destructor
184 : virtual ~SkyEquation();
185 :
186 : // Assignment operator
187 : SkyEquation& operator=(const SkyEquation& other);
188 :
189 : // Copy constructor
190 : SkyEquation(const SkyEquation& other);
191 :
192 : // Add Jones matrices to be used for Gain calculations. Each SkyJones
193 : // knows what type it is and therefore where to slot in.
194 : void setSkyJones(SkyJones& j);
195 :
196 : // Make an approximate PSF for each model. The PSF is approximate
197 : // in the sense that it is a shift-invariant approximation
198 : virtual void makeApproxPSF(casacore::Int model, casacore::ImageInterface<casacore::Float>& PSF);
199 :
200 : // make all the approx psfs in one go
201 : virtual void makeApproxPSF(casacore::PtrBlock<casacore::ImageInterface<casacore::Float> *>& PSFs);
202 :
203 : // Make complex XFRs needed for incrementGradientChiSquared
204 : virtual void makeComplexXFRs();
205 :
206 : // Predict model coherence for the SkyModel. If this is
207 : // incremental then the model visibilities are not reset
208 : // but are simply added to
209 : //virtual void predict(casacore::Bool incremental=false);
210 : virtual void predict(casacore::Bool incremental=false, casacore::MS::PredefinedColumns Type=casacore::MS::MODEL_DATA);
211 :
212 : // Find sum of weights, Chi-squared, and the first and second derivatives
213 : // by transforming to the measurements.
214 : // <group>
215 : virtual void gradientsChiSquared(const casacore::Matrix<casacore::Bool>& required, SkyJones& sj);
216 : virtual void gradientsChiSquared(casacore::Bool incremental, casacore::Bool commitToMS=false);
217 : // </group>
218 :
219 : // Solve for variables. Both the SkyModel and the SkyJones can in
220 : // principle be solved for.
221 : // <group>
222 : virtual casacore::Bool solveSkyModel();
223 : virtual casacore::Bool solveSkyJones(SkyJones& sj);
224 : // </group>
225 :
226 : // Set image plane weighting
227 14 : void setImagePlaneWeighting(const casacore::String& type, const casacore::Float minPB,
228 : const casacore::Float constPB)
229 14 : {scaleType_p = type; minPB_p = minPB; constPB_p = constPB;}
230 :
231 : // Lock and unlock the underlying MeasurementSet
232 : virtual void lock();
233 : virtual void unlock();
234 :
235 : // Return the name of the underlying MeasurementSet
236 : virtual casacore::String associatedMSName();
237 :
238 :
239 : //assign the flux scale that the ftmachines have if they have
240 : virtual void getCoverageImage(casacore::Int model, casacore::ImageInterface<casacore::Float>& im);
241 : //Set this to true if the residual image for mosaic is to be in
242 : //pb^2 units (optimum mode for clean search for centimetric imaging)
243 14 : virtual void doFlatNoise(casacore::Bool doFlat=false){doflat_p=doFlat;};
244 :
245 : //get the weight image from the ftmachines
246 0 : virtual void getWeightImage(const casacore::Int, casacore::ImageInterface<casacore::Float>&){};
247 :
248 0 : virtual casacore::Bool isNewFTM(){return false;}
249 :
250 : virtual void setPhaseCenterTime(const casacore::Double time);
251 : virtual casacore::Double getPhaseCenterTime();
252 : protected:
253 :
254 : // Increment gradientsChiSquared. The image of SkyModel must contain
255 : // the increment to the image. For each model, a collection of
256 : // complex transfer functions are used to avoid gridding and
257 : // degridding all the visibilities.
258 : virtual void incrementGradientsChiSquared();
259 :
260 : // Do the full calculation - this must be called if the FTMachine
261 : // cannot be represented by a Fourier transform
262 :
263 : virtual void fullGradientsChiSquared(casacore::Bool incremental=false);
264 :
265 : SkyEquation() {};
266 :
267 : // Check for validity
268 : casacore::Bool ok();
269 :
270 : // Apply Sky Jones to an image, also adjoint operation
271 : // <group>
272 : virtual void initializeGet(const VisBuffer& vb, casacore::Int row, casacore::Int model,
273 : casacore::Bool incremental);
274 :
275 :
276 : virtual VisBuffer& get(VisBuffer& vb, casacore::Int model, casacore::Bool incremental, casacore::MS::PredefinedColumns Type=casacore::MS::MODEL_DATA);
277 : virtual void finalizeGet();
278 : virtual void initializePut(const VisBuffer &vb, casacore::Int model);
279 :
280 : virtual void put(const VisBuffer& vb, casacore::Int model, casacore::Bool dopsf=false, FTMachine::Type col=FTMachine::OBSERVED);
281 :
282 : virtual void finalizePut(const VisBuffer& vb, casacore::Int Model);
283 :
284 : // This encapsulates all of the change logic we should have to deal
285 : // with (short of returning a range of rows that has the same
286 : // SkyJones). First we look to see if the first row of the VB
287 : // requires a new SkyJones; then we determine if there are
288 : // internal SkyJones changes in the VB which require going
289 : // row by row in the get/put formalism.
290 : virtual void changedSkyJonesLogic(const VisBuffer& vb,
291 : casacore::Bool& firstOneChanges,
292 : casacore::Bool& internalChanges);
293 :
294 : // Have the SkyJones changed since their last application?
295 : virtual casacore::Bool changedSkyJones(const VisBuffer& vb, casacore::Int row);
296 :
297 : // Has the FTMachine changed since last application?
298 : // <group>
299 : virtual casacore::Bool changedFTMachine(const VisBuffer& vb);
300 : virtual casacore::Bool changedIFTMachine(const VisBuffer& vb);
301 : // </group>
302 :
303 : // Do the Sky Jones change in this Visbuffer, starting from row1?
304 : // Returns row2, the last row with the same skyJones
305 : virtual casacore::Bool changedSkyJonesBuffer(const VisBuffer& vb, casacore::Int row1, casacore::Int& row2);
306 :
307 : virtual void resetSkyJones();
308 : virtual void assertSkyJones(const VisBuffer& vb, casacore::Int row);
309 : virtual casacore::ImageInterface<casacore::Complex>& applySkyJones(const VisBuffer& vb, casacore::Int row,
310 : casacore::ImageInterface<casacore::Float>& in,
311 : casacore::ImageInterface<casacore::Complex>& out);
312 : virtual void applySkyJonesInv(const VisBuffer& vb, casacore::Int row,
313 : casacore::ImageInterface<casacore::Complex>& in,
314 : casacore::ImageInterface<casacore::Float>& work,
315 : casacore::ImageInterface<casacore::Float>& out);
316 : virtual void applySkyJonesSquare(const VisBuffer& vb, casacore::Int row,
317 : casacore::Matrix<casacore::Float>& weights,
318 : casacore::ImageInterface<casacore::Float>& work,
319 : casacore::ImageInterface<casacore::Float>& ggS);
320 : // </group>
321 :
322 : // Puts for calculating the complex XFRs
323 : // <group>
324 : virtual void initializePutXFR(const VisBuffer &vb, casacore::Int model, casacore::Int numXFR);
325 : virtual void putXFR(const VisBuffer& vb, casacore::Int model, casacore::Int& numXFR);
326 : virtual void finalizePutXFR(const VisBuffer& vb, casacore::Int Model, casacore::Int numXFR);
327 : // </group>
328 :
329 : // Puts for calculating the complex convolutions
330 : // <group>
331 : virtual void initializePutConvolve(const VisBuffer &vb, casacore::Int model,
332 : casacore::Int numXFR);
333 : virtual void putConvolve(const VisBuffer& vb, casacore::Int model, casacore::Int& numXFR);
334 : virtual void finalizePutConvolve(const VisBuffer& vb, casacore::Int Model, casacore::Int numXFR);
335 : // </group>
336 :
337 : // Get, etc. for a SkyComponent is much simpler
338 : // <group>
339 : virtual VisBuffer& get(VisBuffer& vb, const SkyComponent& component);
340 : virtual VisBuffer& get(VisBuffer& vb, const ComponentList& components);
341 : // Do the sum of the gets for all the models for this visbuffer
342 :
343 : SkyComponent& applySkyJones(SkyComponent& corruptedComponent,
344 : const VisBuffer& vb,
345 : casacore::Int row);
346 : // </group>
347 :
348 : // Modify the ggS and Create the imageScale
349 : virtual void fixImageScale();
350 :
351 : // Deal with scaling or unscaling image or deltaImage
352 : // <group>
353 : virtual void scaleImage(casacore::Int model, casacore::Bool incremental);
354 : virtual void unScaleImage(casacore::Int model, casacore::Bool incremental);
355 : virtual void scaleImage(casacore::Int model=0);
356 : virtual void unScaleImage(casacore::Int model=0);
357 : virtual void scaleDeltaImage(casacore::Int model=0);
358 : virtual void unScaleDeltaImage(casacore::Int model=0);
359 : // </group>
360 :
361 : // Check the VisIter chunck size...force a more reasonable chunk
362 : // if default is too small
363 :
364 : virtual void checkVisIterNumRows(ROVisibilityIterator& vi);
365 :
366 : //virtual void predictComponents(casacore::Bool& incremental, casacore::Bool& initialized);
367 : virtual void predictComponents(casacore::Bool& incremental, casacore::Bool& initialized, casacore::MS::PredefinedColumns Type=casacore::MS::MODEL_DATA);
368 :
369 :
370 : // SkyModel
371 : SkyModel* sm_;
372 :
373 : // VisSet
374 : VisSet* vs_;
375 : //Visibilityiterators
376 : VisibilityIterator* wvi_p;
377 : ROVisibilityIterator* rvi_p;
378 :
379 :
380 : // FTMachine objects
381 : // <group>
382 : FTMachine* ft_;
383 : FTMachine* ift_;
384 : ComponentFTMachine* cft_;
385 : // </group>
386 :
387 : // casacore::List of terms in left to right order
388 : // <group>
389 : SkyJones* ej_;
390 : SkyJones* dj_;
391 : SkyJones* tj_;
392 : SkyJones* fj_;
393 : // </group>
394 :
395 : // Workspace
396 : // <group>
397 : casacore::Float chisq, sumwt;
398 : casacore::Float ggSMax_p;
399 : // </group>
400 :
401 : casacore::LogSink logSink_p;
402 : casacore::LogSink& logSink() {return logSink_p;};
403 :
404 : casacore::Int iDebug_p;
405 :
406 : // Previous VisBuffer, used to determine how to apply
407 : // SkyJones;
408 : // Set in initializePut and initializePutXFR,
409 : // Used in finalizePut and finalizePutXFR
410 :
411 : mutable VisBufferAutoPtr vb_p;
412 :
413 : casacore::Float minPB_p; // ignore model flux below this level in the generalized PB
414 : casacore::Float constPB_p; // make the fluxscale constant for PB above this level
415 : casacore::String scaleType_p; // types: NONE, or SAULT
416 :
417 : casacore::Bool isPSFWork_p; // working for PSF estimation
418 :
419 : casacore::Bool noModelCol_p;
420 :
421 : casacore::Vector<casacore::Bool> modelIsEmpty_p;
422 :
423 : //SkyJones::changed returns a true the first time its called
424 : //We have to ignore this at the very begining and first call to 'changed'
425 : //and not call finalizePut
426 : casacore::Bool isBeginingOfSkyJonesCache_p;
427 : casacore::Bool doflat_p;
428 : };
429 :
430 : } //# NAMESPACE CASA - END
431 :
432 : #endif
433 :
434 :
435 :
436 :
437 :
|
