Line data Source code
1 : //# DOdeconvolver.cc: this implements the deconvolver DO
2 : //# Copyright (C) 1996,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 : //#
27 : //# $Id: DOdeconvolver.cc,v 19.16 2005/12/06 20:18:50 wyoung Exp $
28 :
29 : #include <casacore/casa/Arrays/Matrix.h>
30 : #include <casacore/casa/Arrays/ArrayMath.h>
31 :
32 : #include <casacore/casa/Logging.h>
33 : #include <casacore/casa/Logging/LogIO.h>
34 : #include <casacore/casa/OS/File.h>
35 : #include <casacore/casa/Containers/Record.h>
36 :
37 : #include <casacore/tables/TaQL/TableParse.h>
38 : #include <casacore/tables/Tables/TableRecord.h>
39 : #include <casacore/tables/Tables/TableDesc.h>
40 : #include <casacore/tables/Tables/TableLock.h>
41 : #include <casacore/tables/TaQL/ExprNode.h>
42 :
43 : #include <casacore/casa/BasicSL/String.h>
44 : #include <casacore/casa/Utilities/Assert.h>
45 : #include <casacore/casa/Utilities/Fallible.h>
46 :
47 : #include <casacore/casa/BasicSL/Constants.h>
48 :
49 : #include <casacore/casa/Logging/LogSink.h>
50 : #include <casacore/casa/Logging/LogMessage.h>
51 :
52 : #include <casacore/casa/Arrays/ArrayMath.h>
53 :
54 : #include <casacore/measures/Measures/Stokes.h>
55 : #include <casacore/casa/Quanta/UnitMap.h>
56 : #include <casacore/casa/Quanta/UnitVal.h>
57 : #include <casacore/casa/Quanta/MVAngle.h>
58 : #include <components/ComponentModels/ComponentList.h>
59 : #include <components/ComponentModels/GaussianShape.h>
60 : #include <casacore/measures/Measures/MDirection.h>
61 : #include <casacore/measures/Measures/MPosition.h>
62 : #include <casacore/casa/Quanta/MVEpoch.h>
63 : #include <casacore/measures/Measures/MEpoch.h>
64 :
65 : #include <synthesis/TransformMachines/StokesImageUtil.h>
66 : #include <casacore/lattices/LEL/LatticeExpr.h>
67 : #include <casacore/lattices/LatticeMath/LatticeFFT.h>
68 : #include <casacore/lattices/LatticeMath/LatticeCleaner.h>
69 : #include <casacore/lattices/LatticeMath/LatticeCleanProgress.h>
70 : #include <casacore/lattices/LatticeMath/LatticeConvolver.h>
71 : #include <casacore/lattices/Lattices/TiledLineStepper.h>
72 : #include <casacore/lattices/Lattices/LatticeStepper.h>
73 : #include <casacore/lattices/Lattices/LatticeNavigator.h>
74 : #include <casacore/lattices/Lattices/LatticeIterator.h>
75 : #include <casacore/lattices/Lattices/SubLattice.h>
76 : #include <casacore/lattices/LRegions/LCBox.h>
77 : #include <casacore/lattices/LRegions/LCSlicer.h>
78 :
79 : #include <imageanalysis/ImageAnalysis/ComponentImager.h>
80 : #include <casacore/images/Images/TempImage.h>
81 : #include <casacore/images/Images/PagedImage.h>
82 : #include <casacore/images/Images/ImageSummary.h>
83 : #include <casacore/images/Images/SubImage.h>
84 : #include <casacore/images/Regions/ImageRegion.h>
85 : #include <casacore/images/Images/ImageRegrid.h>
86 : #include <casacore/images/Images/ImageInfo.h>
87 :
88 : #include <casacore/coordinates/Coordinates/CoordinateSystem.h>
89 : #include <casacore/coordinates/Coordinates/DirectionCoordinate.h>
90 : #include <casacore/coordinates/Coordinates/SpectralCoordinate.h>
91 : #include <casacore/coordinates/Coordinates/StokesCoordinate.h>
92 : #include <casacore/coordinates/Coordinates/Projection.h>
93 :
94 : #include <synthesis/MeasurementComponents/ClarkCleanImageSkyModel.h>
95 : #include <synthesis/MeasurementEquations/CEMemProgress.h>
96 : #include <synthesis/MeasurementEquations/CEMemModel.h>
97 : #include <synthesis/MeasurementEquations/ClarkCleanLatModel.h>
98 : #include <synthesis/MeasurementEquations/ClarkCleanProgress.h>
99 : #include <synthesis/DataSampling/ImageDataSampling.h>
100 : #include <synthesis/MeasurementEquations/LatConvEquation.h>
101 : #include <synthesis/MeasurementEquations/IPLatConvEquation.h>
102 : #include <synthesis/MeasurementEquations/Deconvolver.h>
103 : #include <synthesis/MeasurementEquations/Imager.h>
104 : #include <synthesis/MeasurementEquations/ImageMSCleaner.h>
105 : #include <synthesis/MeasurementEquations/ImageNACleaner.h>
106 :
107 :
108 : #include <sstream>
109 :
110 : #include <casacore/casa/namespace.h>
111 : // Implementation comment:
112 : // There are two different philosophies active here:
113 : // ClarkCleanLatModel and CEMemModel are evolutionarily related and
114 : // "solve" a LinearEquation (in this case, a LatticeConvolutionEqaution).
115 : // The LatticeCleaner, which performs Hogbom and MultiScale Cleans,
116 : // has no knowledge of the LatticeConvolutionEqaution (LatConvEquation),
117 : // but carries the PSF and DIRTY around inside itself.
118 : using namespace casacore;
119 : namespace casa { //# NAMESPACE CASA - BEGIN
120 :
121 0 : Deconvolver::Deconvolver()
122 0 : : dirty_p(0),
123 0 : psf_p(0),
124 0 : convolver_p(0),
125 0 : cleaner_p( ),
126 0 : mt_nterms_p(-1),
127 0 : mt_cleaner_p(),
128 0 : mt_valid_p(false)
129 : {
130 :
131 0 : defaults();
132 0 : };
133 :
134 0 : void Deconvolver::defaults()
135 : {
136 0 : mode_p="none";
137 0 : beamValid_p=false;
138 0 : scalesValid_p=false;
139 0 : beam_p = GaussianBeam();
140 0 : residEqn_p = 0;
141 0 : latConvEqn_p = 0;
142 0 : cleaner_p = 0;
143 0 : dirtyName_p = "";
144 0 : psfName_p = "";
145 0 : nx_p=0; ny_p=0; npol_p=0; nchan_p=0;
146 0 : fullPlane_p=false;
147 0 : }
148 :
149 0 : Deconvolver::Deconvolver(const String& dirty, const String& psf)
150 0 : : dirty_p(0), psf_p(0), convolver_p(0), cleaner_p( ), naCleaner_p(nullptr),
151 0 : mt_nterms_p(-1), mt_cleaner_p(), mt_valid_p(false)
152 : {
153 0 : LogIO os(LogOrigin("Deconvolver", "Deconvolver(String& dirty, Strong& psf)", WHERE));
154 0 : defaults();
155 0 : open(dirty, psf);
156 0 : }
157 :
158 0 : Deconvolver::Deconvolver(const Deconvolver &other)
159 0 : : dirty_p(0), psf_p(0), convolver_p(0), cleaner_p( ), naCleaner_p(nullptr),
160 0 : mt_nterms_p(-1), mt_cleaner_p(), mt_valid_p(false)
161 : {
162 0 : defaults();
163 0 : open(other.dirty_p->table().tableName(), other.psf_p->table().tableName());
164 0 : }
165 :
166 0 : Deconvolver &Deconvolver::operator=(const Deconvolver &other)
167 : {
168 0 : if(this != &other){
169 0 : if (dirty_p ) {
170 0 : *dirty_p = *(other.dirty_p);
171 : }
172 0 : if (psf_p ) {
173 0 : *psf_p = *(other.psf_p);
174 : }
175 0 : if (convolver_p && this != &other) {
176 0 : *convolver_p = *(other.convolver_p);
177 : }
178 0 : if ((!cleaner_p.null()) ) {
179 0 : *cleaner_p = *(other.cleaner_p);
180 : }
181 0 : if(naCleaner_p)
182 0 : *naCleaner_p=*(other.naCleaner_p);
183 : }
184 0 : return *this;
185 : }
186 :
187 0 : Deconvolver::~Deconvolver()
188 : {
189 0 : if (psf_p) {
190 0 : delete psf_p;
191 : }
192 0 : psf_p = 0;
193 0 : if (convolver_p) {
194 0 : delete convolver_p;
195 : }
196 0 : convolver_p = 0;
197 : /*if (cleaner_p) {
198 : delete cleaner_p;
199 : }
200 : cleaner_p = 0;
201 : */
202 0 : if (dirty_p) {
203 0 : delete dirty_p;
204 : }
205 0 : dirty_p = 0;
206 0 : }
207 :
208 0 : Bool Deconvolver::open(const String& dirty, const String& psf, Bool warn)
209 : {
210 0 : LogIO os(LogOrigin("Deconvolver", "open()", WHERE));
211 :
212 :
213 : try {
214 0 : if (dirty_p) delete dirty_p; dirty_p = 0;
215 0 : dirty_p = new PagedImage<Float>(dirty);
216 0 : AlwaysAssert(dirty_p, AipsError);
217 0 : nx_p=dirty_p->shape()(0);
218 0 : ny_p=dirty_p->shape()(1);
219 0 : dirty_p->setMaximumCacheSize(2*nx_p*ny_p);
220 0 : dirty_p->setCacheSizeInTiles(10000);
221 :
222 0 : if(dirty_p->shape().nelements()==3){
223 0 : findAxes();
224 0 : if (chanAxis_p > 0)
225 0 : nchan_p=dirty_p->shape()(chanAxis_p);
226 : else
227 0 : nchan_p=0;
228 : }
229 0 : if(dirty_p->shape().nelements()==4){
230 0 : findAxes();
231 0 : npol_p=dirty_p->shape()(polAxis_p);
232 0 : nchan_p=dirty_p->shape()(chanAxis_p);
233 : }
234 0 : dirtyName_p = dirty_p->table().tableName();
235 :
236 0 : if (psf_p) delete psf_p; psf_p = 0;
237 0 : if (psf == ""){
238 0 : if(warn) {
239 : os << LogIO::WARN
240 0 : << "No psf given; please define one before deconvolving" << LogIO::POST;
241 : os << LogIO::WARN
242 0 : << "Use the function open with the psf" << LogIO::POST;
243 : }
244 0 : return true;
245 : }
246 : else{
247 0 : psf_p = new PagedImage<Float>(psf);
248 0 : AlwaysAssert(psf_p, AipsError);
249 0 : psfName_p = psf_p->table().tableName();
250 0 : psf_p->setMaximumCacheSize(2*nx_p*ny_p);
251 0 : psf_p->setCacheSizeInTiles(10000);
252 :
253 : try{
254 0 : os << "Fitting PSF" << LogIO::POST;
255 0 : fitpsf(psf, beam_p);
256 0 : if(! beam_p.isNull()) {
257 0 : os << " Fitted beam is valid"<< LogIO::POST;
258 : }
259 : else {
260 : os << LogIO::WARN << "Fitted beam is invalid: please set using setbeam"
261 0 : << LogIO::POST;
262 : }
263 0 : beamValid_p=true;
264 :
265 0 : } catch (AipsError x) {
266 : os << LogIO::WARN << "Fitted beam is invalid: please set using setbeam"
267 0 : << LogIO::POST;
268 0 : }
269 :
270 0 : if((psf_p->shape()(0) != nx_p) || psf_p->shape()(1) != ny_p){
271 :
272 : os << LogIO::SEVERE
273 0 : << "PSF and Image does not have the same XY shape" << LogIO::POST;
274 : os << LogIO::SEVERE
275 : << "You may wish to regrid the PSF to the same shape as the dirty image"
276 0 : << LogIO::POST;
277 0 : return false;
278 :
279 : }
280 :
281 : try {
282 0 : os << "Making Lattice convolver" << LogIO::POST;
283 0 : if (convolver_p) {
284 0 : delete convolver_p;
285 : }
286 :
287 : // convolver_p = new LatticeConvolver<Float>(*psf_p);
288 : // AlwaysAssert(convolver_p, AipsError);
289 :
290 0 : if (residEqn_p) {
291 0 : delete residEqn_p;
292 : }
293 0 : residEqn_p = 0;
294 :
295 0 : if (latConvEqn_p) {
296 0 : delete latConvEqn_p;
297 : }
298 0 : latConvEqn_p = 0;
299 :
300 0 : os << "Making Image cleaner" << LogIO::POST;
301 : //if (cleaner_p) delete cleaner_p;
302 0 : cleaner_p= new ImageMSCleaner(*psf_p, *dirty_p);
303 0 : naCleaner_p=make_shared<ImageNACleaner>(*psf_p, *dirty_p);
304 0 : if(nchan_p<=1){
305 0 : convolver_p = new LatticeConvolver<Float>(*psf_p);
306 : }
307 : else{
308 0 : if(npol_p > 0 ){
309 0 : IPosition blc(4, 0, 0, 0, 0);
310 0 : IPosition trc(4, nx_p-1, ny_p-1, 0, 0);
311 0 : trc(polAxis_p)=npol_p-1;
312 0 : Slicer sl(blc, trc, Slicer::endIsLast);
313 0 : SubImage<Float> psfSub(*psf_p, sl, true);
314 0 : convolver_p = new LatticeConvolver<Float>(psfSub);
315 0 : AlwaysAssert(convolver_p, AipsError);
316 0 : }
317 : else{
318 0 : IPosition blc(3, 0, 0, 0);
319 0 : IPosition trc(3, nx_p-1, ny_p-1, 0);
320 0 : Slicer sl(blc, trc, Slicer::endIsLast);
321 : //SubImage<Float> dirtySub(*dirty_p, sl, true);
322 0 : SubImage<Float> psfSub(*psf_p, sl, true);
323 0 : convolver_p = new LatticeConvolver<Float>(psfSub);
324 0 : AlwaysAssert(convolver_p, AipsError);
325 0 : }
326 :
327 : }
328 0 : AlwaysAssert(!cleaner_p.null(), AipsError);
329 :
330 0 : return true;
331 :
332 0 : } catch (AipsError x) {
333 0 : os << LogIO::SEVERE << "Caught Exception: "<< x.getMesg() << LogIO::POST;
334 0 : return false;
335 0 : }
336 : }
337 : }
338 0 : catch (AipsError y){
339 0 : throw(y);
340 0 : }
341 0 : }
342 :
343 0 : Bool Deconvolver::reopen()
344 : {
345 0 : LogIO os(LogOrigin("Deconvolver", "reopen()", WHERE));
346 : try {
347 0 : if (dirtyName_p != "" && psfName_p != "") {
348 0 : return (open(dirtyName_p, psfName_p));
349 : } else {
350 0 : return false;
351 : }
352 0 : } catch (AipsError x) {
353 0 : dirty_p->table().unlock();
354 0 : psf_p->table().unlock();
355 0 : os << LogIO::SEVERE << "Caught Exception: "<< x.getMesg() << LogIO::POST;
356 0 : return false;
357 0 : }
358 : return false;
359 0 : }
360 :
361 : // Fit the psf. If psf is blank then make the psf first.
362 0 : Bool Deconvolver::fitpsf(const String& psf, GaussianBeam& beam)
363 : {
364 0 : if(!valid()) return false;
365 :
366 0 : LogIO os(LogOrigin("Deconvolver", "fitpsf()", WHERE));
367 :
368 : try {
369 :
370 0 : os << "Fitting to psf" << LogIO::POST;
371 :
372 0 : if(psf=="") {
373 0 : os << LogIO::SEVERE << "Need a psf name" << LogIO::POST;
374 0 : return false;
375 : }
376 :
377 0 : PagedImage<Float> psfImage(psf);
378 0 : StokesImageUtil::FitGaussianPSF(psfImage, beam);
379 0 : beam_p = beam;
380 0 : beamValid_p=true;
381 :
382 : os << " Beam fit: " << beam_p.getMajor("arcsec") << " by "
383 : << beam_p.getMinor("arcsec") << " (arcsec) at pa "
384 0 : << beam_p.getPA(Unit("deg")) << " (deg) " << endl;
385 0 : return true;
386 0 : } catch (AipsError x) {
387 0 : os << LogIO::SEVERE << "Exception: " << x.getMesg() << LogIO::POST;
388 0 : }
389 0 : return true;
390 0 : }
391 :
392 0 : Bool Deconvolver::close()
393 : {
394 0 : if(!valid()) return false;
395 0 : if (detached()) return true;
396 0 : LogIO os(LogOrigin("Deconvolver", "close()", WHERE));
397 :
398 0 : os << "Closing images and detaching from Deconvolver" << LogIO::POST;
399 0 : if(psf_p) delete psf_p; psf_p = 0;
400 0 : if(dirty_p) delete dirty_p; dirty_p = 0;
401 0 : if (convolver_p) delete convolver_p; convolver_p = 0;
402 0 : if (residEqn_p) delete residEqn_p; residEqn_p = 0;
403 0 : if (latConvEqn_p) delete latConvEqn_p; latConvEqn_p = 0;
404 :
405 0 : return true;
406 0 : }
407 :
408 0 : String Deconvolver::dirtyname() const
409 : {
410 0 : if (detached()) {
411 0 : return "none";
412 : }
413 0 : return dirty_p->table().tableName();
414 : }
415 :
416 0 : String Deconvolver::psfname() const
417 : {
418 0 : if (detached()) {
419 0 : return "none";
420 : }
421 0 : return psf_p->table().tableName();
422 : }
423 :
424 0 : Bool Deconvolver::summary() const
425 : {
426 0 : if(!valid()) return false;
427 0 : LogOrigin OR("Deconvolver", "Deconvolver::summary()", WHERE);
428 :
429 0 : LogIO los(OR);
430 :
431 : try {
432 :
433 0 : los << "Summary of dirty image" << LogIO::POST;
434 0 : dirty_p->table().lock();
435 : {
436 0 : ImageSummary<Float> ims(*dirty_p);
437 0 : ims.list(los);
438 0 : }
439 :
440 0 : los << endl << state() << LogIO::POST;
441 0 : dirty_p->table().unlock();
442 :
443 0 : los << "Summary of PSF" << LogIO::POST;
444 0 : psf_p->table().lock();
445 : {
446 0 : ImageSummary<Float> psfs(*psf_p);
447 0 : psfs.list(los);
448 0 : }
449 :
450 0 : los << "Summary of scales" << LogIO::POST;
451 0 : if(scalesValid_p) {
452 0 : los << "Scales set" << LogIO::POST;
453 : }
454 : else {
455 0 : los << "Scales not set" << LogIO::POST;
456 : }
457 :
458 0 : los << endl << state() << LogIO::POST;
459 0 : psf_p->table().unlock();
460 0 : return true;
461 0 : } catch (AipsError x) {
462 : los << LogIO::SEVERE << "Caught Exception: " << x.getMesg()
463 0 : << LogIO::POST;
464 0 : dirty_p->table().unlock();
465 0 : psf_p->table().unlock();
466 0 : return false;
467 0 : }
468 :
469 : return true;
470 0 : }
471 :
472 0 : String Deconvolver::state() const
473 : {
474 0 : ostringstream os;
475 :
476 : try {
477 0 : os << "General: " << endl;
478 0 : if(dirty_p != 0){
479 0 : os << " Dirty image is " << dirty_p->table().tableName() << endl;
480 0 : dirty_p->table().unlock();
481 : }
482 0 : if(psf_p !=0){
483 0 : os << " PSF is " << psf_p->table().tableName() << endl;
484 0 : psf_p->table().unlock();
485 : }
486 0 : if(beamValid_p) {
487 0 : os << " Beam fit: " << beam_p.getMajor("arcsec") << " by "
488 0 : << beam_p.getMinor("arcsec") << " (arcsec) at pa "
489 0 : << beam_p.getPA(Unit("deg")) << " (deg) " << endl;
490 : }
491 : else {
492 0 : os << " Beam fit is not valid" << endl;
493 : }
494 :
495 0 : } catch (AipsError x) {
496 0 : LogOrigin OR("Deconvolver", "Deconvolver::state()", WHERE);
497 0 : LogIO los(OR);
498 : los << LogIO::SEVERE << "Caught exception: " << x.getMesg()
499 0 : << LogIO::POST;
500 0 : dirty_p->table().unlock();
501 0 : psf_p->table().unlock();
502 0 : }
503 0 : return String(os);
504 0 : }
505 :
506 : // Restore: at least one model must be supplied
507 0 : Bool Deconvolver::restore(const String& model, const String& image,
508 : GaussianBeam& mbeam)
509 : {
510 :
511 0 : if(!valid()) return false;
512 0 : LogIO os(LogOrigin("Deconvolver", "restore()", WHERE));
513 :
514 0 : dirty_p->table().lock();
515 0 : psf_p->table().lock();
516 : try {
517 :
518 : // Validate the names
519 0 : if(model=="") {
520 : os << LogIO::SEVERE << "Need a model"
521 0 : << LogIO::POST;
522 0 : return false;
523 : }
524 :
525 0 : String imagename(image);
526 0 : if(imagename=="") imagename=model+".restored";
527 0 : removeTable(imagename);
528 0 : if(!clone(model, imagename)) return false;
529 :
530 : // Smooth all the images and add residuals
531 0 : PagedImage<Float> modelImage0(model);
532 :
533 : //
534 0 : TiledShape tShape(dirty_p->shape());
535 0 : ImageInterface<Float>* modelImage_p = new TempImage<Float>(tShape, dirty_p->coordinates());
536 : //
537 0 : ImageRegrid<Float> regridder;
538 0 : Vector<Double> locate;
539 0 : Bool missedIt = regridder.insert(*modelImage_p, locate, modelImage0);
540 : //cerr << "missedIt " << missedIt << endl;
541 0 : if (!missedIt) {
542 0 : os << LogIO::SEVERE << "Problem in getting model Image on correct grid " << LogIO::POST;
543 : }
544 :
545 0 : PagedImage<Float> imageImage(modelImage_p->shape(),
546 : modelImage_p->coordinates(),
547 0 : image);
548 :
549 0 : TempImage<Float> dirtyModelImage(modelImage_p->shape(),modelImage_p->coordinates());
550 : //imageImage.copyData(*modelImage_p);
551 0 : if(! mbeam.isNull()) {
552 : os << " Using specified beam: " << mbeam.getMajor("arcsec") << " by "
553 : << mbeam.getMinor("arcsec") << " (arcsec) at pa "
554 0 : << mbeam.getPA(Unit("deg")) << " (deg) " << LogIO::POST;
555 : //StokesImageUtil::Convolve(imageImage, mbeam, false);
556 : }
557 : else {
558 0 : if(! beam_p.isNull()) {
559 : os << " Using fitted beam: " << beam_p.getMajor("arcsec") << " by "
560 : << beam_p.getMinor("arcsec") << " (arcsec) at pa "
561 0 : << beam_p.getPA(Unit("deg")) << " (deg) " << LogIO::POST;
562 : //StokesImageUtil::Convolve(imageImage, beam_p, false);
563 0 : mbeam = beam_p;
564 : }
565 : else {
566 0 : os << LogIO::SEVERE << "Restoring beam not specified" << LogIO::POST;
567 0 : return false;
568 : }
569 : }
570 : //Model * restoring beam
571 : {
572 0 : IPosition convshp=modelImage_p->shape();
573 0 : convshp[0]=nx_p; convshp[1]=ny_p;
574 0 : for (uInt k=2; k< convshp.nelements(); ++k) convshp[k]=1;
575 0 : TempImage<Float> gaussim(convshp, modelImage_p->coordinates());
576 0 : gaussim.set(0.0);
577 0 : ImageInfo ii = gaussim.imageInfo();
578 0 : ii.setRestoringBeam(mbeam);
579 0 : gaussim.setImageInfo(ii);
580 0 : gaussim.setUnits(Unit("Jy/beam"));
581 0 : putGaussian(gaussim, mbeam);
582 : //////////////////////
583 : //PagedImage<Float>xx(gaussim.shape(), gaussim.coordinates(), "tempGauss");
584 : //xx.copyData(gaussim);
585 : //////////////////
586 0 : LatticeConvolver<Float> lc(gaussim);
587 0 : lc.linear(imageImage, *modelImage_p);
588 0 : }
589 : // PSF * Model
590 0 : convolver_p->circular(dirtyModelImage, *modelImage_p);
591 :
592 : // Smoothed + Dirty - PSF * Model
593 0 : imageImage.copyData(LatticeExpr<Float>(imageImage+*dirty_p - dirtyModelImage));
594 : {
595 0 : ImageInfo ii = imageImage.imageInfo();
596 0 : ii.setRestoringBeam(mbeam);
597 0 : imageImage.setImageInfo(ii);
598 0 : imageImage.setUnits(Unit("Jy/beam"));
599 0 : }
600 :
601 0 : dirty_p->table().unlock();
602 0 : psf_p->table().unlock();
603 0 : if (modelImage_p != & modelImage0) {
604 0 : delete modelImage_p;
605 : }
606 0 : return true;
607 0 : } catch (AipsError x) {
608 0 : dirty_p->table().unlock();
609 : os << LogIO::SEVERE << "Caught exception: " << x.getMesg()
610 0 : << LogIO::POST;
611 0 : }
612 0 : dirty_p->table().unlock();
613 0 : return true;
614 0 : }
615 :
616 : // Residual
617 0 : Bool Deconvolver::residual(const String& model, const String& image)
618 : {
619 :
620 0 : if(!valid()) return false;
621 0 : LogIO os(LogOrigin("Deconvolver", "residual()", WHERE));
622 :
623 0 : dirty_p->table().lock();
624 0 : psf_p->table().lock();
625 : try {
626 :
627 : // Validate the names
628 0 : if(model=="") {
629 : os << LogIO::SEVERE << "Need a model"
630 0 : << LogIO::POST;
631 : }
632 :
633 0 : String imagename(image);
634 0 : if(imagename=="") imagename=model+".residual";
635 0 : removeTable(imagename);
636 0 : if(!clone(dirty_p->table().tableName(), imagename)) return false;
637 :
638 : // Smooth all the images and add residuals
639 :
640 : // modelImage_p is a pointer to an image with the model data in it, but the
641 : // shape of the dirty image
642 0 : PagedImage<Float> modelImage0(model);
643 :
644 0 : TiledShape tShape(dirty_p->shape());
645 0 : ImageInterface<Float>* modelImage_p = new TempImage<Float>(tShape, dirty_p->coordinates());
646 : //
647 0 : ImageRegrid<Float> regridder;
648 0 : Vector<Double> locate;
649 0 : Bool missedIt = regridder.insert(*modelImage_p, locate, modelImage0);
650 0 : if (!missedIt) {
651 0 : os << LogIO::SEVERE << "Problem in getting model Image on correct grid " << LogIO::POST;
652 : }
653 :
654 0 : PagedImage<Float> imageImage(modelImage_p->shape(),
655 : modelImage_p->coordinates(),
656 0 : image);
657 : // PSF * Model
658 0 : convolver_p->circular(imageImage, *modelImage_p);
659 :
660 : // Dirty - PSF * Model
661 0 : imageImage.copyData(LatticeExpr<Float>(*dirty_p-imageImage));
662 0 : imageImage.setUnits(Unit("Jy/beam"));
663 :
664 0 : dirty_p->table().unlock();
665 0 : psf_p->table().unlock();
666 0 : if (modelImage_p != & modelImage0) {
667 0 : delete modelImage_p;
668 : }
669 0 : return true;
670 0 : } catch (AipsError x) {
671 0 : dirty_p->table().unlock();
672 0 : psf_p->table().unlock();
673 : os << LogIO::SEVERE << "Caught exception: " << x.getMesg()
674 0 : << LogIO::POST;
675 0 : }
676 0 : dirty_p->table().unlock();
677 0 : psf_p->table().unlock();
678 0 : return true;
679 0 : }
680 :
681 : // Make an empty image
682 0 : Bool Deconvolver::make(const String& model)
683 : {
684 :
685 0 : if(!valid()) return false;
686 0 : LogIO os(LogOrigin("Deconvolver", "make()", WHERE));
687 :
688 0 : dirty_p->table().lock();
689 : try {
690 :
691 : // Make an image with the required shape and coordinates
692 0 : String modelName(model);
693 0 : if(modelName=="") modelName=dirty_p->table().tableName()+".model";
694 0 : os << "Making empty image: " << model << LogIO::POST;
695 :
696 0 : removeTable(modelName);
697 0 : PagedImage<Float> modelImage(dirty_p->shape(),
698 0 : dirty_p->coordinates(), model);
699 0 : modelImage.set(0.0);
700 :
701 0 : modelImage.table().tableInfo().setSubType("GENERIC");
702 0 : modelImage.setUnits(Unit("Jy/pixel"));
703 0 : dirty_p->table().unlock();
704 0 : return true;
705 0 : } catch (AipsError x) {
706 0 : dirty_p->table().unlock();
707 0 : os << LogIO::SEVERE << "Exception: " << x.getMesg() << LogIO::POST;
708 0 : }
709 0 : dirty_p->table().unlock();
710 0 : return true;
711 0 : };
712 :
713 :
714 : // Make an empty image, but with only ONE STOKES pixel
715 0 : Bool Deconvolver::make1(const String& model)
716 : {
717 :
718 0 : if(!valid()) return false;
719 0 : LogIO os(LogOrigin("Deconvolver", "make1()", WHERE));
720 :
721 0 : dirty_p->table().lock();
722 : try {
723 :
724 : // Make an image with the required shape and coordinates
725 0 : String modelName(model);
726 0 : if(modelName=="") modelName=dirty_p->table().tableName()+".model";
727 0 : os << "Making empty image: " << model << LogIO::POST;
728 :
729 0 : removeTable(modelName);
730 0 : IPosition newshape = dirty_p->shape();
731 0 : newshape(2) = 1;
732 : PagedImage<Float> modelImage(newshape,
733 0 : dirty_p->coordinates(), model);
734 0 : modelImage.set(0.0);
735 :
736 0 : modelImage.table().tableInfo().setSubType("GENERIC");
737 0 : modelImage.setUnits(Unit("Jy/pixel"));
738 0 : dirty_p->table().unlock();
739 0 : return true;
740 0 : } catch (AipsError x) {
741 0 : dirty_p->table().unlock();
742 0 : os << LogIO::SEVERE << "Exception: " << x.getMesg() << LogIO::POST;
743 0 : }
744 0 : dirty_p->table().unlock();
745 0 : return true;
746 0 : };
747 :
748 :
749 : // Make an empty image, modeled in templateImage
750 0 : Bool Deconvolver::make(const String& model, ImageInterface<Float>& templateImage)
751 : {
752 :
753 0 : if(!valid()) return false;
754 0 : LogIO os(LogOrigin("Deconvolver", "make()", WHERE));
755 :
756 : try {
757 :
758 : // Make an image with the required shape and coordinates
759 0 : String modelName(model);
760 :
761 0 : os << "Making empty image: " << model << LogIO::POST;
762 0 : removeTable(modelName);
763 0 : PagedImage<Float> modelImage(templateImage.shape(),
764 0 : templateImage.coordinates(), model);
765 0 : modelImage.set(0.0);
766 :
767 0 : modelImage.table().tableInfo().setSubType("GENERIC");
768 0 : modelImage.setUnits(Unit("Jy/pixel"));
769 0 : return true;
770 0 : } catch (AipsError x) {
771 0 : os << LogIO::SEVERE << "Exception: " << x.getMesg() << LogIO::POST;
772 0 : }
773 0 : return true;
774 0 : };
775 :
776 :
777 : SubImage<Float>*
778 0 : Deconvolver::innerQuarter(PagedImage<Float>& in)
779 : {
780 :
781 0 : IPosition blc(in.shape().nelements(), 0);
782 0 : IPosition trc(in.shape()-1);
783 0 : for (Int i=0;i<Int(in.shape().nelements());i++) {
784 0 : blc(i)=in.shape()(i)/4;
785 0 : trc(i)=blc(i)+in.shape()(i)/2-1;
786 0 : if(trc(i)<0) trc(i)=1;
787 : }
788 0 : LCSlicer quarter(blc, trc);
789 0 : SubImage<Float>* si = new SubImage<Float>(in, quarter, true);
790 0 : return si;
791 0 : };
792 :
793 :
794 : SubImage<Float>*
795 0 : Deconvolver::allQuarters(PagedImage<Float>& in)
796 : {
797 0 : SubImage<Float>* si = new SubImage<Float>(in, true);
798 0 : return si;
799 : };
800 :
801 0 : Bool Deconvolver::smooth(const String& model,
802 : const String& image,
803 : GaussianBeam& mbeam,
804 : Bool normalizeVolume)
805 : {
806 0 : if(!valid()) return false;
807 0 : LogIO os(LogOrigin("Deconvolver", "smooth()", WHERE));
808 :
809 0 : dirty_p->table().lock();
810 : try {
811 :
812 0 : os << "Smoothing image" << LogIO::POST;
813 :
814 0 : if(model=="") {
815 0 : os << LogIO::SEVERE << "Need a name for model " << LogIO::POST;
816 0 : return false;
817 : }
818 :
819 0 : if(mbeam.getMajor().getValue()==0.0) {
820 0 : if(beamValid_p) {
821 0 : os << "Using previous beam fit" << LogIO::POST;
822 0 : mbeam = beam_p;
823 :
824 : }
825 : else {
826 0 : os << LogIO::SEVERE << "Specified beam is invalid" << LogIO::POST;
827 : }
828 : }
829 :
830 : // Smooth all the images
831 0 : PagedImage<Float> modelImage(model);
832 0 : PagedImage<Float> imageImage(modelImage.shape(),
833 : modelImage.coordinates(),
834 0 : image);
835 : //
836 0 : imageImage.copyData(modelImage);
837 0 : StokesImageUtil::Convolve(imageImage, mbeam, normalizeVolume);
838 : {
839 0 : ImageInfo ii = imageImage.imageInfo();
840 0 : ii.setRestoringBeam(mbeam);
841 0 : imageImage.setImageInfo(ii);
842 0 : imageImage.setUnits(Unit("Jy/beam"));
843 0 : }
844 :
845 0 : dirty_p->table().unlock();
846 0 : psf_p->table().unlock();
847 0 : return true;
848 0 : } catch (AipsError x) {
849 0 : dirty_p->table().unlock();
850 0 : psf_p->table().unlock();
851 0 : os << LogIO::SEVERE << "Exception: " << x.getMesg() << LogIO::POST;
852 0 : }
853 0 : dirty_p->table().unlock();
854 0 : psf_p->table().unlock();
855 0 : return true;
856 0 : }
857 : //will clean casa axes order images
858 0 : Bool Deconvolver::clarkclean(const Int niter,
859 : const Float gain, const Quantity& threshold,
860 : const String& model, const String& maskName,
861 : Float& maxresidual, Int& iterused,
862 : Float cycleFactor){
863 :
864 0 : Bool retval=false;
865 0 : Double thresh=threshold.get("Jy").getValue();
866 0 : String imagename(model);
867 : // Make first image with the required shape and coordinates only if
868 : // it doesn't exist yet. Otherwise we'll throw an exception later
869 0 : if(imagename=="") imagename=dirty_p->table().tableName()+".clarkclean";
870 0 : if(!Table::isWritable(imagename)) {
871 0 : make(imagename);
872 : }
873 0 : PagedImage<Float> modelImage(imagename);
874 0 : Bool hasMask=false;
875 0 : if(maskName !="")
876 0 : hasMask=Table::isReadable(maskName);
877 0 : if(hasMask){
878 0 : PagedImage<Float> mask(maskName);
879 0 : retval=ClarkCleanImageSkyModel::clean(modelImage, *dirty_p, *psf_p, mask, maxresidual, iterused, gain, niter, thresh, cycleFactor, true, true);
880 0 : }
881 : else{
882 0 : ImageInterface<Float> *tmpMask=0;
883 0 : retval=ClarkCleanImageSkyModel::clean(modelImage, *dirty_p, *psf_p, *tmpMask, maxresidual, iterused, gain, niter, thresh, cycleFactor, false, true);
884 : }
885 :
886 0 : return retval;
887 0 : }
888 : // Clean algorithm
889 0 : Bool Deconvolver::clarkclean(const Int niter,
890 : const Float gain, const Quantity& threshold,
891 : const Bool displayProgress,
892 : const String& model, const String& maskName,
893 : const Int histBins,
894 : const Vector<Int>& vi_psfPatchSize, const Float maxExtPsf,
895 : const Float speedUp, Int maxNumPix,
896 : const Int maxNumMajorCycles,
897 : const Int maxNumMinorIterations)
898 : {
899 0 : if(!valid()) return false;
900 0 : LogIO os(LogOrigin("Deconvolver", "clarkclean()", WHERE));
901 :
902 0 : IPosition psfPatchSize(vi_psfPatchSize);
903 :
904 :
905 0 : dirty_p->table().lock();
906 0 : psf_p->table().lock();
907 :
908 0 : String imagename(model);
909 : // Make first image with the required shape and coordinates only if
910 : // it doesn't exist yet. Otherwise we'll throw an exception later
911 0 : if(imagename=="") imagename=dirty_p->table().tableName()+".clarkclean";
912 0 : if(!Table::isWritable(imagename)) {
913 0 : make(imagename);
914 : }
915 0 : PagedImage<Float> modelImage(imagename);
916 0 : ClarkCleanProgress *ccpp = 0;
917 :
918 : {
919 0 : ostringstream oos;
920 0 : oos << "Clean gain = " <<gain<<", Niter = "<<niter<<", Threshold = "
921 0 : <<threshold; ;
922 0 : os << String(oos) << LogIO::POST;
923 0 : }
924 : {
925 0 : ostringstream oos;
926 0 : oos << "nHhistBins = "
927 0 : <<histBins << ", maxExtPsf = "<<maxExtPsf<<", psfPatchSize = "
928 0 : <<psfPatchSize<<", maxNumPix = "<<maxNumPix;
929 0 : os << String(oos) << LogIO::POST;
930 0 : }
931 : {
932 0 : ostringstream oos;
933 0 : oos << "Speedup Factor = "<<speedUp<<", maxMajorCycles = "
934 0 : << maxNumMajorCycles<<", maxMinorIterations = "<<maxNumMinorIterations;
935 0 : os << String(oos) << LogIO::POST;
936 0 : }
937 :
938 0 : os << "Cleaning image using Clark Clean algorithm" << LogIO::POST;
939 :
940 :
941 :
942 :
943 :
944 : try {
945 0 : if(model=="") {
946 0 : os << LogIO::SEVERE << "Need a name for model " << LogIO::POST;
947 0 : return false;
948 : }
949 0 : PagedImage<Float> *mask = 0;
950 0 : Bool isCubeMask=false;
951 :
952 : Int xbeg, xend, ybeg, yend;
953 : //default clean box
954 0 : xbeg=nx_p/4;
955 0 : xend=3*nx_p/4-1;
956 0 : ybeg=ny_p/4;
957 0 : yend=3*ny_p/4-1;
958 :
959 : // Deal with mask
960 0 : if (maskName != "") {
961 0 : if( Table::isReadable(maskName)) {
962 0 : mask= new PagedImage<Float>(maskName);
963 0 : if (chanAxis_p < Int(mask->shape().nelements())){
964 0 : if (mask->shape()(chanAxis_p) > 1)
965 0 : isCubeMask=true;
966 : }
967 0 : checkMask(*mask, xbeg, xend, ybeg, yend);
968 0 : AlwaysAssert(mask, AipsError);
969 :
970 : } else {
971 0 : os << LogIO::SEVERE << "Mask "<< mask<<" is not readable" << LogIO::POST;
972 : }
973 : }
974 0 : SubImage<Float> *maskSub=0;
975 0 : IPosition blc(2,xbeg,ybeg);
976 0 : IPosition trc(2,xend,yend);
977 0 : Bool result=false;
978 0 : if(nchan_p >=1){
979 0 : for (Int k=0; k<nchan_p; ++k){
980 0 : os<< "Cleaning channel " << k+1 << LogIO::POST;
981 0 : if(npol_p > 0 ){
982 0 : blc.resize(4);
983 0 : blc(chanAxis_p)=k;
984 0 : blc(polAxis_p)=0;
985 0 : trc.resize(4);
986 0 : trc(polAxis_p)=npol_p-1;
987 0 : trc(chanAxis_p)=k;
988 :
989 : }
990 : else{
991 0 : blc.resize(3);
992 0 : trc.resize(3);
993 0 : blc(chanAxis_p)=k;
994 0 : trc(chanAxis_p)=k;
995 :
996 : }
997 :
998 :
999 :
1000 0 : Slicer sl(blc, trc, Slicer::endIsLast);
1001 0 : SubImage<Float> psfSub;
1002 :
1003 0 : if(mask != 0){
1004 0 : if( (isCubeMask) || (!isCubeMask && maskSub == 0 )){
1005 0 : if(maskSub !=0 ) delete maskSub;
1006 0 : blc(0)=0; blc(1)=0;
1007 0 : trc(0)=nx_p-1; trc(1)=ny_p-1;
1008 0 : sl=Slicer(blc, trc, Slicer::endIsLast);
1009 0 : maskSub=new SubImage<Float> (*mask,sl,false);
1010 0 : checkMask(*maskSub, xbeg, xend, ybeg, yend);
1011 0 : blc(0)=xbeg; blc(1)=ybeg;
1012 0 : trc(0)=xend; trc(1)=yend;
1013 0 : sl =Slicer(blc, trc, Slicer::endIsLast);
1014 0 : delete maskSub;
1015 0 : maskSub=new SubImage<Float> (*mask,sl,false);
1016 : }
1017 : }
1018 :
1019 :
1020 0 : SubImage<Float> dirtySub(*dirty_p, sl, true);
1021 0 : SubImage<Float> modelSub(modelImage,sl,true);
1022 0 : IPosition blc_psf=blc; IPosition trc_psf=trc;
1023 0 : if(psf_p->shape().nelements() != dirty_p->shape().nelements()){
1024 0 : blc_psf.resize(psf_p->shape().nelements());
1025 0 : trc_psf.resize(psf_p->shape().nelements());
1026 : }
1027 0 : blc_psf(0)=0; blc_psf(1)=0;
1028 0 : trc_psf(0)=nx_p-1; trc_psf(1)=ny_p-1;
1029 0 : sl=Slicer(blc_psf, trc_psf, Slicer::endIsLast);
1030 0 : psfSub=SubImage<Float>(*psf_p, sl, true);
1031 :
1032 0 : ClarkCleanLatModel myClarkCleaner(modelSub);
1033 0 : if(mask !=0 )
1034 0 : myClarkCleaner.setMask(*maskSub);
1035 :
1036 0 : myClarkCleaner.setNumberIterations(niter);
1037 0 : if (maxNumMajorCycles > 0 )
1038 0 : myClarkCleaner.setMaxNumberMajorCycles((uInt)maxNumMajorCycles);
1039 0 : if (maxNumMinorIterations > 0 )
1040 0 : myClarkCleaner.setMaxNumberMinorIterations((uInt)maxNumMinorIterations);
1041 :
1042 0 : myClarkCleaner.setGain(gain);
1043 0 : Double d_thresh = threshold.getValue("Jy");
1044 0 : myClarkCleaner.setThreshold((Float)d_thresh);
1045 :
1046 0 : myClarkCleaner.setPsfPatchSize(psfPatchSize);
1047 0 : myClarkCleaner.setHistLength((uInt)histBins);
1048 0 : myClarkCleaner.setMaxExtPsf(maxExtPsf);
1049 0 : myClarkCleaner.setSpeedup(speedUp);
1050 :
1051 0 : if (maxNumPix == 0)
1052 0 : maxNumPix = (Int)(modelImage.shape().product()*0.04);
1053 0 : myClarkCleaner.setMaxNumPix((uInt)maxNumPix);
1054 :
1055 : //Now actually do the clean
1056 0 : if (displayProgress) {
1057 0 : ccpp = new ClarkCleanProgress ();
1058 0 : myClarkCleaner.setProgress(*ccpp);
1059 : }
1060 0 : if(latConvEqn_p !=0) delete latConvEqn_p;
1061 0 : latConvEqn_p=0;
1062 0 : latConvEqn_p = new LatConvEquation (psfSub, dirtySub);
1063 0 : result=myClarkCleaner.solve(*latConvEqn_p);
1064 0 : }
1065 : }
1066 : else{
1067 0 : IPosition blc(modelImage.shape().nelements(),0);
1068 0 : Int elem= npol_p >0 ? npol_p:0;
1069 0 : IPosition trc(modelImage.shape().nelements(),elem);
1070 0 : blc(0)=xbeg; blc(1)=ybeg;
1071 0 : trc(0)=xend; trc(1)=yend;
1072 0 : Slicer sl(blc, trc, Slicer::endIsLast);
1073 0 : SubImage<Float> maskSub;
1074 0 : SubImage<Float> dirtySub(*dirty_p, sl, true);
1075 0 : SubImage<Float> modelSub(modelImage,sl,true);
1076 0 : if(psf_p->shape().nelements() != dirty_p->shape().nelements()){
1077 0 : blc.resize(psf_p->shape().nelements());
1078 0 : trc.resize(psf_p->shape().nelements());
1079 : }
1080 0 : blc(0)=0; blc(1)=0;
1081 0 : trc(0)=nx_p-1; trc(1)=ny_p-1;
1082 0 : sl=Slicer(blc, trc, Slicer::endIsLast);
1083 0 : SubImage<Float> psfSub(*psf_p, sl, true);
1084 :
1085 0 : ClarkCleanLatModel myClarkCleaner(modelSub);
1086 0 : if(mask !=0 ){
1087 0 : maskSub= SubImage<Float>(*mask, sl, false);
1088 0 : myClarkCleaner.setMask(maskSub);
1089 : }
1090 :
1091 0 : myClarkCleaner.setNumberIterations(niter);
1092 0 : if (maxNumMajorCycles > 0 )
1093 0 : myClarkCleaner.setMaxNumberMajorCycles((uInt)maxNumMajorCycles);
1094 0 : if (maxNumMinorIterations > 0 )
1095 0 : myClarkCleaner.setMaxNumberMinorIterations((uInt)maxNumMinorIterations);
1096 :
1097 0 : myClarkCleaner.setGain(gain);
1098 0 : Double d_thresh = threshold.getValue("Jy");
1099 0 : myClarkCleaner.setThreshold((Float)d_thresh);
1100 :
1101 0 : myClarkCleaner.setPsfPatchSize(psfPatchSize);
1102 0 : myClarkCleaner.setHistLength((uInt)histBins);
1103 0 : myClarkCleaner.setMaxExtPsf(maxExtPsf);
1104 0 : myClarkCleaner.setSpeedup(speedUp);
1105 :
1106 0 : if (maxNumPix == 0)
1107 0 : maxNumPix = (Int)(modelImage.shape().product()*0.04);
1108 0 : myClarkCleaner.setMaxNumPix((uInt)maxNumPix);
1109 :
1110 : //Now actually do the clean
1111 0 : if (displayProgress) {
1112 0 : ccpp = new ClarkCleanProgress ();
1113 0 : myClarkCleaner.setProgress(*ccpp);
1114 : }
1115 0 : latConvEqn_p = new LatConvEquation (psfSub, dirtySub);
1116 0 : result=myClarkCleaner.solve(*latConvEqn_p);
1117 0 : }
1118 0 : dirty_p->table().unlock();
1119 0 : psf_p->table().unlock();
1120 0 : if (ccpp != 0) { delete ccpp; ccpp = 0;}
1121 0 : delete latConvEqn_p; latConvEqn_p = 0;
1122 0 : if (mask) { delete mask;}
1123 :
1124 0 : return result;
1125 0 : } catch (AipsError x) {
1126 0 : os << LogIO::SEVERE << "Exception: " << x.getMesg() << LogIO::POST;
1127 0 : }
1128 0 : dirty_p->table().unlock();
1129 0 : psf_p->table().unlock();
1130 0 : if (ccpp != 0) {delete ccpp; ccpp = 0; }
1131 0 : delete latConvEqn_p; latConvEqn_p = 0;
1132 :
1133 :
1134 0 : return true;
1135 0 : };
1136 :
1137 :
1138 :
1139 0 : Bool Deconvolver::setupLatCleaner(const String& /*algorithm*/, const Int /*niter*/,
1140 : const Float /*gain*/, const Quantity& /*threshold*/,
1141 : const Bool /*displayProgress*/){
1142 :
1143 0 : LogIO os(LogOrigin("Deconvolver", "clean()", WHERE));
1144 :
1145 : /*
1146 : if((algorithm=="msclean")||(algorithm=="fullmsclean" || algorithm=="multiscale" || algorithm=="fullmultiscale")) {
1147 : os << "Cleaning image using multi-scale algorithm" << LogIO::POST;
1148 : if(!scalesValid_p) {
1149 : os << LogIO::SEVERE << "Scales not yet set" << LogIO::POST;
1150 : return false;
1151 : }
1152 : cleaner_p->setscales(scaleSizes_p);
1153 : cleaner_p->setcontrol(CleanEnums::MULTISCALE, niter, gain, threshold);
1154 : }
1155 : else if (algorithm=="hogbom") {
1156 : if(!scalesValid_p) {
1157 : Vector<Float> dummy;
1158 : setscales("nscales", 1, dummy);
1159 : }
1160 : cleaner_p->setscales(scaleSizes_p);
1161 : cleaner_p->setcontrol(CleanEnums::HOGBOM, niter, gain, threshold);
1162 : } else {
1163 : os << LogIO::SEVERE << "Unknown algorithm: " << algorithm << LogIO::POST;
1164 : return false;
1165 : }
1166 :
1167 :
1168 :
1169 :
1170 : if(algorithm=="fullmsclean" || algorithm=="fullmultiscale") {
1171 : os << "Cleaning full image using multi-scale algorithm" << LogIO::POST;
1172 : cleaner_p->ignoreCenterBox(true);
1173 : }
1174 : */
1175 0 : return true;
1176 :
1177 0 : }
1178 : // Clean algorithm
1179 0 : Bool Deconvolver::clean(const String& algorithm, const Int niter,
1180 : const Float gain, const Quantity& threshold,
1181 : const Bool displayProgress,
1182 : const String& model, const String& mask, Float& maxResidual, Int& iterationsDone)
1183 : {
1184 :
1185 0 : if(!valid()) return false;
1186 0 : LogIO os(LogOrigin("Deconvolver", "clean()", WHERE));
1187 :
1188 0 : dirty_p->table().lock();
1189 0 : psf_p->table().lock();
1190 : try {
1191 :
1192 0 : if(model=="") {
1193 0 : os << LogIO::SEVERE << "Need a name for model " << LogIO::POST;
1194 0 : return false;
1195 : }
1196 : //Int psfnchan=psf_p->shape()(chanAxis_p);
1197 : //Int masknchan=0;
1198 :
1199 0 : String imagename(model);
1200 : // Make first image with the required shape and coordinates only if
1201 : // it doesn't exist yet. Otherwise we'll throw an exception later
1202 0 : if(imagename=="") imagename=dirty_p->table().tableName()+"."+algorithm;
1203 0 : if(!Table::isWritable(imagename)) {
1204 0 : make(imagename);
1205 : }
1206 :
1207 : {
1208 0 : ostringstream oos;
1209 0 : oos << "Clean gain = " <<gain<<", Niter = "<<niter<<", Threshold = "
1210 0 : <<threshold << ", Algorithm " << algorithm ;
1211 0 : os << String(oos) << LogIO::POST;
1212 0 : }
1213 :
1214 0 : PagedImage<Float> modelImage(imagename);
1215 :
1216 0 : AlwaysAssert(!cleaner_p.null(), AipsError);
1217 0 : PagedImage<Float> *maskim = 0;
1218 : // Deal with mask
1219 0 : if (mask != "") {
1220 0 : if( Table::isReadable(mask)) {
1221 0 : maskim = new PagedImage<Float>(mask);
1222 0 : AlwaysAssert(maskim, AipsError);
1223 0 : cleaner_p->setMask(*maskim);
1224 : } else {
1225 :
1226 0 : os << LogIO::SEVERE << "Mask "<< mask<<" is not readable" << LogIO::POST;
1227 : }
1228 : }
1229 :
1230 :
1231 0 : Bool result=false;
1232 :
1233 0 : result=cleaner_p->clean(modelImage, algorithm, niter, gain, threshold, displayProgress);
1234 0 : maxResidual=cleaner_p->maxResidual();
1235 0 : iterationsDone=cleaner_p->numberIterations();
1236 0 : dirty_p->table().relinquishAutoLocks(true);
1237 0 : dirty_p->table().unlock();
1238 0 : psf_p->table().relinquishAutoLocks(true);
1239 0 : psf_p->table().unlock();
1240 0 : if (maskim) delete maskim;
1241 :
1242 0 : return result;
1243 0 : } catch (AipsError x) {
1244 0 : dirty_p->table().unlock();
1245 0 : psf_p->table().unlock();
1246 0 : os << LogIO::SEVERE << "Exception: " << x.getMesg() << LogIO::POST;
1247 0 : }
1248 :
1249 0 : return true;
1250 0 : }
1251 :
1252 0 : Bool Deconvolver::naclean(const Int niter,
1253 : const Float gain, const Quantity& threshold,
1254 : const String& model, const String& maskname, const Int masksupp, const Int memType, const Float numSigma, Float& maxResidual, Int& iterationsDone)
1255 : {
1256 :
1257 0 : if(!valid()) return false;
1258 0 : LogIO os(LogOrigin("Deconvolver", "clean()", WHERE));
1259 :
1260 0 : dirty_p->table().lock();
1261 0 : psf_p->table().lock();
1262 : try {
1263 :
1264 0 : if(model=="") {
1265 0 : os << LogIO::SEVERE << "Need a name for model " << LogIO::POST;
1266 0 : return false;
1267 : }
1268 : //Int psfnchan=psf_p->shape()(chanAxis_p);
1269 : //Int masknchan=0;
1270 :
1271 0 : String imagename(model);
1272 : // Make first image with the required shape and coordinates only if
1273 : // it doesn't exist yet. Otherwise we'll throw an exception later
1274 0 : if(imagename=="") imagename=dirty_p->table().tableName()+".naclean";
1275 0 : if(!Table::isWritable(imagename)) {
1276 0 : make(imagename);
1277 : }
1278 :
1279 : {
1280 0 : ostringstream oos;
1281 0 : oos << "Clean gain = " <<gain<<", Niter = "<<niter<<", Threshold = "
1282 0 : <<threshold ;
1283 0 : os << String(oos) << LogIO::POST;
1284 0 : }
1285 :
1286 0 : PagedImage<Float> modelImage(imagename);
1287 :
1288 0 : AlwaysAssert(naCleaner_p != nullptr, AipsError);
1289 0 : PagedImage<Float> *maskim = 0;
1290 : // Deal with mask
1291 0 : String mask=maskname;
1292 0 : if (mask == "") {
1293 0 : mask=dirty_p->table().tableName()+".mask";
1294 : }
1295 0 : if(!Table::isWritable(mask))
1296 0 : make(mask);
1297 :
1298 0 : maskim = new PagedImage<Float>(mask);
1299 :
1300 :
1301 0 : AlwaysAssert(maskim, AipsError);
1302 0 : naCleaner_p->setMask(*maskim);
1303 :
1304 :
1305 :
1306 :
1307 0 : Bool result=false;
1308 :
1309 0 : result=naCleaner_p->clean(modelImage, niter, gain, threshold, masksupp, memType, numSigma);
1310 0 : maxResidual=cleaner_p->maxResidual();
1311 0 : iterationsDone=cleaner_p->iteration();
1312 0 : dirty_p->table().relinquishAutoLocks(true);
1313 0 : dirty_p->table().unlock();
1314 0 : psf_p->table().relinquishAutoLocks(true);
1315 0 : psf_p->table().unlock();
1316 0 : if (maskim) delete maskim;
1317 :
1318 0 : return result;
1319 0 : } catch (AipsError x) {
1320 0 : dirty_p->table().unlock();
1321 0 : psf_p->table().unlock();
1322 0 : os << LogIO::SEVERE << "Exception: " << x.getMesg() << LogIO::POST;
1323 0 : }
1324 :
1325 0 : return true;
1326 0 : }
1327 :
1328 :
1329 :
1330 : // MEM algorithm
1331 0 : Bool Deconvolver::mem(const String& entropy, const Int niter,
1332 : const Quantity& sigma, const Quantity& targetFlux,
1333 : Bool constrainTargetFlux, Bool displayProgress,
1334 : const String& model, const String& priorImage,
1335 : const String& maskImage,
1336 : const Bool imagePlane)
1337 : {
1338 :
1339 0 : if(!valid()) return false;
1340 0 : LogIO os(LogOrigin("Deconvolver", "mem()", WHERE));
1341 :
1342 : // SubImage<Float>* dirtyQ_p =0;
1343 : // SubImage<Float>* modelImageQ_p =0;
1344 : // SubImage<Float>* priorImageQ_p =0;
1345 : // SubImage<Float>* maskImageQ_p =0;
1346 :
1347 0 : Entropy* myEnt_p =0;
1348 0 : CEMemProgress * memProgress_p = 0;
1349 0 : residEqn_p=0;
1350 :
1351 : try {
1352 :
1353 : // if(model=="") {
1354 : // os << LogIO::SEVERE << "Need a name for model " << LogIO::POST;
1355 : // return false;
1356 : // }
1357 :
1358 0 : Bool initializeModel = false;
1359 : Int xbeg, xend, ybeg, yend;
1360 0 : if (imagePlane) {
1361 0 : xbeg=0;
1362 0 : xend=nx_p-1;
1363 0 : ybeg=0;
1364 0 : yend=ny_p-1;
1365 0 : fullPlane_p=true;
1366 : } else {
1367 0 : xbeg=nx_p/4;
1368 0 : xend=3*nx_p/4-1;
1369 0 : ybeg=ny_p/4;
1370 0 : yend=3*ny_p/4-1;
1371 : }
1372 :
1373 0 : String imagename(model);
1374 : // Make first image with the required shape and coordinates only if
1375 : // it doesn't exist yet. Otherwise we'll throw an exception later
1376 0 : if(imagename=="") {
1377 0 : imagename=dirty_p->table().tableName()+"."+entropy;
1378 : os << LogIO::WARN << "No model name given, model will be "
1379 0 : << imagename << LogIO::POST;
1380 : }
1381 0 : if(!Table::isWritable(imagename)) {
1382 0 : initializeModel = true;
1383 0 : make(imagename);
1384 0 : dirty_p->table().lock();
1385 0 : psf_p->table().lock();
1386 : }
1387 :
1388 0 : PagedImage<Float> modelImage(imagename);
1389 :
1390 : // if (imagePlane) {
1391 : // modelImageQ_p = allQuarters(modelImage);
1392 : // } else {
1393 : // modelImageQ_p = innerQuarter(modelImage);
1394 : // }
1395 :
1396 : {
1397 0 : ostringstream oos;
1398 0 : oos << "MEM Niter = "<<niter<<", Sigma = "<<sigma <<
1399 0 : ", TargetFlux = " <<targetFlux <<
1400 0 : ", ConstrainTargetFlux = " <<constrainTargetFlux <<
1401 0 : ", Entropy = " << entropy;
1402 0 : os << String(oos) << LogIO::POST;
1403 0 : }
1404 :
1405 0 : if(entropy=="entropy") {
1406 0 : os << "Deconvolving image using maximum entropy algorithm" << LogIO::POST;
1407 0 : myEnt_p = new EntropyI;
1408 : }
1409 0 : else if (entropy=="emptiness") {
1410 0 : myEnt_p = new EntropyEmptiness;
1411 : }
1412 : else {
1413 0 : os << " Known MEM entropies: entropy | emptiness " << LogIO::POST;
1414 0 : os << LogIO::SEVERE << "Unknown MEM entropy: " << entropy << LogIO::POST;
1415 0 : return false;
1416 : }
1417 :
1418 :
1419 0 : PagedImage<Float> *mask = 0;
1420 0 : Bool isCubeMask=false;
1421 0 : PagedImage<Float> *prior =0;
1422 :
1423 : // Deal with mask
1424 0 : if (maskImage != "") {
1425 0 : if( Table::isReadable(maskImage)) {
1426 0 : mask= new PagedImage<Float>(maskImage);
1427 0 : if (chanAxis_p < Int(mask->shape().nelements())){
1428 0 : if (mask->shape()(chanAxis_p) > 1)
1429 0 : isCubeMask=true;
1430 : }
1431 0 : checkMask(*mask, xbeg, xend, ybeg, yend);
1432 0 : AlwaysAssert(mask, AipsError);
1433 :
1434 : } else {
1435 0 : os << LogIO::SEVERE << "Mask "<< mask<<" is not readable" << LogIO::POST;
1436 : }
1437 : }
1438 :
1439 0 : if (priorImage!="") {
1440 0 : if( Table::isReadable(priorImage)) {
1441 0 : prior= new PagedImage<Float>(priorImage);
1442 : } else {
1443 0 : os << LogIO::SEVERE << "Prior "<< prior<<" is not readable" << LogIO::POST;
1444 : }
1445 : }
1446 :
1447 0 : Bool result=false;
1448 0 : SubImage<Float> *maskSub=0;
1449 0 : IPosition blc(2,xbeg,ybeg);
1450 0 : IPosition trc(2,xend,yend);
1451 0 : if(nchan_p >= 1){
1452 0 : for (Int k=0; k< nchan_p; ++k){
1453 0 : os<< "Processing channel " << k+1 << LogIO::POST;
1454 0 : if(npol_p > 0 ){
1455 0 : blc.resize(4);
1456 0 : blc(chanAxis_p)=k;
1457 0 : blc(polAxis_p)=0;
1458 0 : trc.resize(4);
1459 0 : trc(polAxis_p)=npol_p-1;
1460 0 : trc(chanAxis_p)=k;
1461 :
1462 : }
1463 : else{
1464 0 : blc.resize(3);
1465 0 : trc.resize(3);
1466 0 : blc(chanAxis_p)=k;
1467 0 : trc(chanAxis_p)=k;
1468 :
1469 : }
1470 :
1471 0 : Slicer sl(blc, trc, Slicer::endIsLast);
1472 0 : SubImage<Float> psfSub;
1473 0 : SubImage<Float> priorSub;
1474 :
1475 0 : if(mask != 0){
1476 0 : if( (isCubeMask) || (!isCubeMask && maskSub == 0 )){
1477 0 : if(maskSub !=0 ) delete maskSub;
1478 0 : blc(0)=0; blc(1)=0;
1479 0 : trc(0)=nx_p-1; trc(1)=ny_p-1;
1480 0 : sl=Slicer(blc, trc, Slicer::endIsLast);
1481 0 : maskSub=new SubImage<Float> (*mask,sl,false);
1482 0 : checkMask(*maskSub, xbeg, xend, ybeg, yend);
1483 0 : blc(0)=xbeg; blc(1)=ybeg;
1484 0 : trc(0)=xend; trc(1)=yend;
1485 0 : sl =Slicer(blc, trc, Slicer::endIsLast);
1486 0 : delete maskSub;
1487 0 : maskSub=new SubImage<Float> (*mask,sl,false);
1488 : }
1489 : }
1490 :
1491 :
1492 0 : if(prior !=0 ){
1493 0 : priorSub= SubImage<Float>(*prior, false);
1494 : }
1495 :
1496 0 : SubImage<Float> dirtySub(*dirty_p, sl, true);
1497 0 : SubImage<Float> modelSub(modelImage,sl,true);
1498 0 : IPosition blc_psf=blc; IPosition trc_psf=trc;
1499 0 : if(psf_p->shape().nelements() != dirty_p->shape().nelements()){
1500 0 : blc_psf.resize(psf_p->shape().nelements());
1501 0 : trc_psf.resize(psf_p->shape().nelements());
1502 : }
1503 0 : blc_psf(0)=0; blc_psf(1)=0;
1504 0 : trc_psf(0)=nx_p-1; trc_psf(1)=ny_p-1;
1505 0 : sl=Slicer(blc_psf, trc_psf, Slicer::endIsLast);
1506 0 : psfSub=SubImage<Float>(*psf_p, sl, true);
1507 :
1508 :
1509 :
1510 0 : CEMemModel myMemer( *myEnt_p, modelSub, niter, sigma.getValue("Jy"),
1511 0 : targetFlux.getValue("Jy"), constrainTargetFlux,
1512 0 : initializeModel, imagePlane);
1513 :
1514 0 : if (!initializeModel) {
1515 0 : Record info=modelImage.miscInfo();
1516 : try {
1517 0 : Float alpha = 0.0;
1518 0 : Float beta = 0.0;
1519 0 : info.get("ALPHA", alpha);
1520 0 : myMemer.setAlpha(alpha);
1521 0 : info.get("BETA", beta);
1522 0 : myMemer.setBeta(beta);
1523 0 : } catch (AipsError x) {
1524 : // could not get Alpha and Beta for initialization
1525 : // continue
1526 : os << "Could not retrieve Alpha and Beta from previously initialized model"
1527 0 : << LogIO::POST;
1528 0 : }
1529 0 : }
1530 :
1531 :
1532 :
1533 0 : if(prior != 0){
1534 0 : myMemer.setPrior(priorSub);
1535 : }
1536 0 : if (mask != 0) {
1537 0 : myMemer.setMask(*maskSub);
1538 : }
1539 :
1540 :
1541 : // Now actually do the MEM deconvolution
1542 0 : if (displayProgress) {
1543 0 : memProgress_p = new CEMemProgress ();
1544 0 : myMemer.setProgress(*memProgress_p);
1545 : }
1546 :
1547 0 : if(residEqn_p !=0) delete residEqn_p;
1548 0 : residEqn_p=0;
1549 0 : if (imagePlane) {
1550 0 : residEqn_p = new IPLatConvEquation (psfSub, dirtySub);
1551 : } else {
1552 0 : residEqn_p = new LatConvEquation (psfSub, dirtySub);
1553 : }
1554 :
1555 0 : result=myMemer.solve(*residEqn_p);
1556 :
1557 0 : Record info=modelImage.miscInfo();
1558 0 : info.define("ALPHA", myMemer.getBeta());
1559 0 : info.define("BETA", myMemer.getAlpha());
1560 0 : modelImage.setMiscInfo(info);
1561 0 : }
1562 : }
1563 :
1564 : else{
1565 0 : SubImage<Float>* dirtyQ =0;
1566 0 : SubImage<Float>* modelQ =0;
1567 0 : Bool initializeModel = false;
1568 :
1569 0 : if (imagePlane) {
1570 0 : dirtyQ = allQuarters(*dirty_p);
1571 : } else {
1572 0 : dirtyQ = innerQuarter(*dirty_p);
1573 : }
1574 0 : if (imagePlane) {
1575 0 : modelQ = allQuarters(modelImage);
1576 : } else {
1577 0 : modelQ = innerQuarter(modelImage);
1578 : }
1579 : CEMemModel myMemer( *myEnt_p, *modelQ, niter,
1580 0 : sigma.getValue("Jy"),
1581 0 : targetFlux.getValue("Jy"), constrainTargetFlux,
1582 0 : initializeModel, imagePlane);
1583 0 : if (!initializeModel) {
1584 0 : Record info=modelImage.miscInfo();
1585 : try {
1586 0 : Float alpha = 0.0;
1587 0 : Float beta = 0.0;
1588 0 : info.get("ALPHA", alpha);
1589 0 : myMemer.setAlpha(alpha);
1590 0 : info.get("BETA", beta);
1591 0 : myMemer.setBeta(beta);
1592 0 : } catch (AipsError x) {
1593 : // could not get Alpha and Beta for initialization
1594 : // continue
1595 : os << "Could not retrieve Alpha and Beta from previously initialized model"
1596 0 : << LogIO::POST;
1597 0 : }
1598 0 : }
1599 0 : SubImage<Float> * priorQ=NULL;
1600 0 : if(prior !=0){
1601 :
1602 0 : if (imagePlane) {
1603 0 : priorQ = allQuarters(*prior);
1604 : } else {
1605 0 : priorQ = innerQuarter(*prior);
1606 :
1607 : }
1608 0 : myMemer.setPrior(*priorQ);
1609 : }
1610 0 : SubImage<Float> *maskQ=NULL;
1611 0 : if(mask !=0){
1612 :
1613 0 : if (imagePlane) {
1614 0 : maskQ = allQuarters(*mask);
1615 : } else {
1616 0 : maskQ = innerQuarter(*mask);
1617 :
1618 : }
1619 0 : myMemer.setMask(*maskQ);
1620 : }
1621 0 : if (displayProgress) {
1622 0 : memProgress_p = new CEMemProgress ();
1623 0 : myMemer.setProgress(*memProgress_p);
1624 : }
1625 :
1626 0 : if (imagePlane) {
1627 0 : residEqn_p = new IPLatConvEquation (*psf_p, *dirtyQ);
1628 : } else {
1629 0 : residEqn_p = new LatConvEquation (*psf_p, *dirtyQ);
1630 : }
1631 :
1632 0 : result=myMemer.solve(*residEqn_p);
1633 :
1634 0 : Record info=modelImage.miscInfo();
1635 0 : info.define("ALPHA", myMemer.getBeta());
1636 0 : info.define("BETA", myMemer.getAlpha());
1637 0 : modelImage.setMiscInfo(info);
1638 :
1639 0 : if (dirtyQ != 0) {delete dirtyQ; dirtyQ = 0;}
1640 0 : if (modelQ != 0) {delete modelQ; modelQ = 0;}
1641 0 : if (residEqn_p != 0) {delete residEqn_p; residEqn_p = 0;}
1642 0 : if (priorQ != 0) {delete priorQ; priorQ = 0;}
1643 0 : if (maskQ != 0) {delete maskQ; maskQ = 0;}
1644 :
1645 0 : }
1646 :
1647 0 : modelImage.setUnits(Unit("Jy/pixel"));
1648 :
1649 0 : dirty_p->table().unlock();
1650 0 : psf_p->table().unlock();
1651 0 : if (myEnt_p != 0) {delete myEnt_p; myEnt_p = 0;}
1652 0 : if (memProgress_p!=0) {delete memProgress_p; memProgress_p = 0; }
1653 :
1654 0 : return result;
1655 :
1656 0 : } catch (AipsError x) {
1657 0 : os << LogIO::SEVERE << "Exception: " << x.getMesg() << LogIO::POST;
1658 0 : }
1659 :
1660 0 : dirty_p->table().unlock();
1661 0 : psf_p->table().unlock();
1662 0 : if (myEnt_p != 0) {delete myEnt_p; myEnt_p = 0;}
1663 0 : if (residEqn_p != 0) {delete residEqn_p; residEqn_p = 0;}
1664 0 : if (memProgress_p!=0) {delete memProgress_p; memProgress_p = 0; }
1665 :
1666 0 : return true;
1667 0 : }
1668 :
1669 : // makeprior, for MEM
1670 0 : Bool Deconvolver::makeprior(const String& prior, const String& templatename,
1671 : const Quantity& lowClipFrom, const Quantity& lowClipTo,
1672 : const Quantity& highClipFrom, const Quantity& highClipTo,
1673 : const Vector<Int>& blc, const Vector<Int>& trc)
1674 : {
1675 :
1676 0 : if(!valid()) return false;
1677 0 : LogIO os(LogOrigin("Deconvolver", "makeprior()", WHERE));
1678 :
1679 : try {
1680 0 : if(templatename=="") {
1681 : os << LogIO::SEVERE << "Need a name for template image " << endl
1682 : << "May I suggest you make a clean or mem image and " << endl
1683 0 : << "smooth it for the template?" << LogIO::POST;
1684 0 : return false;
1685 : }
1686 0 : if(prior=="") {
1687 0 : os << LogIO::SEVERE << "Need a name for output prior image " << LogIO::POST;
1688 0 : return false;
1689 : }
1690 :
1691 0 : PagedImage<Float> templateImage(templatename);
1692 0 : String priorname(prior);
1693 0 : if(priorname=="") priorname=templateImage.table().tableName()+".prior";
1694 0 : if(!Table::isWritable(priorname)) {
1695 0 : make(priorname);
1696 : }
1697 : //
1698 0 : PagedImage<Float> priorImage(priorname);
1699 0 : ImageInterface<Float>* templateImage2_p = 0;
1700 : //
1701 0 : if (priorImage.shape() != templateImage.shape()) {
1702 0 : TiledShape tShape(priorImage.shape());
1703 0 : templateImage2_p = new TempImage<Float>(tShape, priorImage.coordinates());
1704 : //
1705 0 : ImageRegrid<Float> regridder;
1706 0 : Vector<Double> locate;
1707 0 : Bool missedIt = regridder.insert(*templateImage2_p, locate, templateImage);
1708 0 : if (!missedIt) {
1709 0 : os << LogIO::SEVERE << "Problem in getting template Image on correct grid " << LogIO::POST;
1710 : }
1711 0 : } else {
1712 0 : templateImage2_p = &templateImage;
1713 : }
1714 :
1715 : {
1716 0 : ostringstream oos;
1717 0 : oos << "Prior = "<<priorname<<", template = "<<templatename << endl;
1718 0 : oos <<" Clip Below = " << lowClipFrom.getValue("Jy") <<
1719 0 : ", Replace with = " << lowClipTo.getValue("Jy") << endl;
1720 :
1721 0 : oos <<" Clip Above = " << highClipFrom.getValue("Jy") <<
1722 0 : ", Replace with = " << highClipTo.getValue("Jy") << endl;
1723 : // oos <<" blc = " << blc <<", trc = " << trc << endl;
1724 0 : os << String(oos) << LogIO::POST;
1725 0 : }
1726 :
1727 0 : priorImage.set(lowClipTo.getValue("Jy"));
1728 :
1729 0 : IPosition iblc(blc);
1730 0 : IPosition itrc(trc);
1731 0 : IPosition imshape(priorImage.shape());
1732 : // Index::convertIPosition(iblc, blc);
1733 : // Index::convertIPosition(itrc, trc);
1734 0 : IPosition iinc(imshape.nelements(),1);
1735 0 : LCBox::verify(iblc, itrc, iinc, imshape);
1736 0 : LCSlicer box(iblc, itrc);
1737 :
1738 0 : SubImage<Float> templateBox(*templateImage2_p, ImageRegion(box), false);
1739 0 : SubImage<Float> priorBox(priorImage, ImageRegion(box), true);
1740 :
1741 : // do Low clipping
1742 0 : priorBox.copyData( (LatticeExpr<Float>)
1743 0 : (iif(templateBox<lowClipFrom.getValue("Jy"),
1744 : lowClipTo.getValue("Jy"), templateBox)) );
1745 : // do High clipping
1746 0 : priorBox.copyData( (LatticeExpr<Float>)
1747 0 : (iif(priorBox > highClipFrom.getValue("Jy"),
1748 : highClipTo.getValue("Jy"), priorBox)) );
1749 :
1750 :
1751 0 : return true;
1752 0 : } catch (AipsError x) {
1753 0 : os << LogIO::SEVERE << "Exception: " << x.getMesg() << LogIO::POST;
1754 0 : }
1755 :
1756 0 : return true;
1757 0 : }
1758 :
1759 : // clipimage
1760 0 : Bool Deconvolver::clipimage(const String& clippedImageName, const String& inputImageName,
1761 : const Quantity& threshold)
1762 : {
1763 0 : if(!valid()) return false;
1764 0 : LogIO os(LogOrigin("Deconvolver", "clipimage()", WHERE));
1765 :
1766 : try {
1767 0 : if(inputImageName=="") {
1768 0 : os << LogIO::SEVERE << "Need a name for the image to clip" << LogIO::POST;
1769 0 : return false;
1770 : }
1771 0 : if(clippedImageName=="") {
1772 0 : os << LogIO::SEVERE << "Need a name for output clipped image " << LogIO::POST;
1773 0 : return false;
1774 : }
1775 :
1776 0 : PagedImage<Float> inputImage(inputImageName);
1777 0 : String clippedImageName2(clippedImageName);
1778 0 : if(clippedImageName2=="") clippedImageName2 = inputImage.table().tableName()+".clipped";
1779 0 : if(!Table::isWritable(clippedImageName2)) {
1780 0 : make (clippedImageName2);
1781 : }
1782 0 : PagedImage<Float> clippedImage(clippedImageName2);
1783 0 : if (clippedImage.shape() != inputImage.shape() ) {
1784 0 : os << LogIO::SEVERE << "Input and clipped image sizes disagree " << LogIO::POST;
1785 0 : return false;
1786 : }
1787 : {
1788 0 : ostringstream oos;
1789 0 : oos << "Clipped Image = "<<clippedImageName2<<", Input Image = "<< inputImageName << endl;
1790 0 : oos << "Clip with Stokes I below = " << threshold.getValue("Jy");
1791 0 : os << String(oos) << LogIO::POST;
1792 0 : }
1793 :
1794 0 : IPosition trc = inputImage.shape() - 1;
1795 0 : IPosition blc(4,0);
1796 :
1797 0 : trc(2) = 0;
1798 0 : blc(2) = 0;
1799 0 : LCSlicer boxI(blc, trc);
1800 0 : SubImage<Float> stokesISub(inputImage, ImageRegion(boxI), false);
1801 : Int iStokes;
1802 0 : for (iStokes=0; iStokes < inputImage.shape()(2); iStokes++) {
1803 0 : trc(2) = iStokes;
1804 0 : blc(2) = iStokes;
1805 0 : LCSlicer box(blc, trc);
1806 0 : SubImage<Float> stokesClippedSub(clippedImage, ImageRegion(box), true);
1807 0 : SubImage<Float> stokesInputSub(inputImage, ImageRegion(box), false);
1808 0 : if (stokesISub.shape() != stokesClippedSub.shape() ) {
1809 0 : os << LogIO::SEVERE << "Input and clipped image sizes disagree " << LogIO::POST;
1810 0 : return false;
1811 : }
1812 0 : stokesClippedSub.copyData( (LatticeExpr<Float>)
1813 0 : (iif(stokesISub < threshold.getValue("Jy"),
1814 : 0.0, stokesInputSub)) );
1815 0 : }
1816 0 : return true;
1817 0 : } catch (AipsError x) {
1818 0 : os << LogIO::SEVERE << "Exception: " << x.getMesg() << LogIO::POST;
1819 0 : }
1820 :
1821 0 : return true;
1822 0 : }
1823 :
1824 : // boxmask
1825 0 : Bool Deconvolver::boxmask(const String& boxmask,
1826 : const Vector<Int> blc, const Vector<Int> trc,
1827 : const Quantity& fillValue, const Quantity& externalValue)
1828 : {
1829 0 : if(!valid()) return false;
1830 0 : LogIO os(LogOrigin("Deconvolver", "boxmask()", WHERE));
1831 :
1832 : try {
1833 0 : if(boxmask=="") {
1834 0 : os << LogIO::SEVERE << "Need a name for output boxmask image " << LogIO::POST;
1835 0 : return false;
1836 : }
1837 0 : String boxname(boxmask);
1838 0 : if(boxname=="") boxname="boxmask";
1839 0 : if(!Table::isWritable(boxname)) {
1840 0 : make(boxname);
1841 : }
1842 0 : PagedImage<Float> boxImage(boxname);
1843 :
1844 : {
1845 0 : ostringstream oos;
1846 0 : oos << "BoxMask = "<<boxname<<
1847 0 : ", blc = " << blc(0) << " " << blc(1)<<
1848 0 : ", trc = " << trc(0) << " " << trc(1);
1849 0 : os << String(oos) << LogIO::POST;
1850 0 : }
1851 :
1852 0 : boxImage.set(externalValue.getValue("Jy"));
1853 :
1854 : // This only makes a 2-d box; will need to fix this for other
1855 : // image sorts
1856 0 : uInt dim = boxImage.ndim();
1857 0 : IPosition pshape = boxImage.shape();
1858 0 : IPosition blc0(dim, 0);
1859 0 : IPosition trc0(dim, 0);
1860 0 : blc0(0) = max(0, blc(0));
1861 0 : blc0(1) = max(0, blc(1));
1862 0 : if (trc0(0) == 0) trc0(0) = pshape(0)-1;
1863 0 : if (trc0(1) == 0) trc0(1) = pshape(1)-1;
1864 0 : trc0(0) = min(trc(0), pshape(0)-1);
1865 0 : trc0(1) = min(trc(1), pshape(1)-1);
1866 0 : LCSlicer box(blc0, trc0);
1867 0 : SubImage<Float> innerSub(boxImage, ImageRegion(box), true);
1868 0 : innerSub.set(fillValue.getValue("Jy"));
1869 0 : return true;
1870 0 : } catch (AipsError x) {
1871 0 : os << LogIO::SEVERE << "Exception: " << x.getMesg() << LogIO::POST;
1872 0 : }
1873 0 : return true;
1874 0 : }
1875 :
1876 : //regionmask
1877 0 : Bool Deconvolver::regionmask(const String& maskimage, Record* imageRegRec, Matrix<Quantity>& blctrcs, const Float& value){
1878 :
1879 0 : LogIO os(LogOrigin("deconvolver", "regionmask()", WHERE));
1880 0 : if(!dirty_p) {
1881 : os << LogIO::SEVERE << "Program logic error: Dirty image pointer dirty_p not yet set"
1882 0 : << LogIO::POST;
1883 0 : return false;
1884 : }
1885 0 : if(!Table::isWritable(maskimage)) {
1886 0 : if (!clone(dirty_p->name(),maskimage)) return false;
1887 0 : PagedImage<Float> mim(maskimage);
1888 0 : mim.set(0.0);
1889 0 : mim.table().relinquishAutoLocks();
1890 0 : }
1891 0 : Matrix<Float> circles;
1892 0 : return Imager::regionToImageMask(maskimage, imageRegRec, blctrcs, circles, value);
1893 :
1894 0 : }
1895 :
1896 : // Fourier transform the model
1897 0 : Bool Deconvolver::ft(const String& model, const String& transform)
1898 : {
1899 0 : if(!valid()) return false;
1900 :
1901 0 : LogIO os(LogOrigin("Deconvolver", "ft()", WHERE));
1902 :
1903 : try {
1904 :
1905 0 : if(model=="") {
1906 0 : os << LogIO::SEVERE << "Need a name for model " << LogIO::POST;
1907 0 : return false;
1908 : }
1909 :
1910 0 : os << "Fourier transforming model" << LogIO::POST;
1911 :
1912 0 : String transformname(transform);
1913 0 : if(transformname=="") transformname=model+".ft";
1914 0 : removeTable(transformname);
1915 :
1916 0 : PagedImage<Float> modelImage(model);
1917 0 : PagedImage<Complex> transformImage(modelImage.shape(),
1918 : modelImage.coordinates(),
1919 0 : transformname);
1920 0 : transformImage.copyData(LatticeExpr<Complex>(toComplex(modelImage)));
1921 :
1922 0 : LatticeFFT::cfft2d(transformImage);
1923 :
1924 0 : return true;
1925 0 : } catch (AipsError x) {
1926 0 : os << LogIO::SEVERE << "Exception: " << x.getMesg() << LogIO::POST;
1927 0 : }
1928 0 : return true;
1929 0 : }
1930 :
1931 0 : Bool Deconvolver::setscales(const String& scaleMethod,
1932 : const Int inscales,
1933 : const Vector<Float>& userScaleSizes)
1934 : {
1935 0 : LogIO os(LogOrigin("Deconvolver", "setscales()", WHERE));
1936 :
1937 0 : AlwaysAssert(!cleaner_p.null(), AipsError);
1938 :
1939 0 : Vector<Double> cells = psf_p->coordinates().increment();
1940 0 : os << "Cell size = " << abs(cells(0)/C::arcsec) << LogIO::POST;
1941 0 : AlwaysAssert (cells(0)!=0.0, AipsError);
1942 0 : scaleSizes_p.resize();
1943 0 : if (scaleMethod == "nscales") {
1944 0 : Int nscales=inscales;
1945 :
1946 0 : if(nscales<1) {
1947 0 : os << "Using default of 5 scales" << LogIO::POST;
1948 0 : nscales=5;
1949 : }
1950 :
1951 : // Validate scales
1952 0 : Float scaleInc=beam_p.getMinor().get("arcsec").getValue()/abs(cells(0)/C::arcsec);
1953 :
1954 0 : Vector<Float> scaleSizes(nscales);
1955 : os << "Creating " << nscales <<
1956 0 : " scales from powerlaw nscales method" << LogIO::POST;
1957 0 : scaleSizes(0) = 0.0;
1958 0 : os << "scale 1 = 0.0 pixels " << LogIO::POST;
1959 0 : for (Int scale=1; scale<nscales;scale++) {
1960 0 : scaleSizes(scale) =
1961 0 : scaleInc * pow(10.0, (Float(scale)-2.0)/2.0);
1962 0 : os << "scale " << scale+1 << " = " << scaleSizes(scale)
1963 0 : << " pixels" << LogIO::POST;
1964 : }
1965 0 : scaleSizes_p=scaleSizes;
1966 0 : cleaner_p->setscales(scaleSizes);
1967 0 : scalesValid_p=true;
1968 :
1969 0 : } else if (scaleMethod == "uservector") {
1970 0 : if (userScaleSizes.nelements() <= 0) {
1971 : os << LogIO::SEVERE
1972 : << "Need at least one scale for method uservector"
1973 0 : << LogIO::POST;
1974 : }
1975 0 : os << "Creating scales from uservector method: " << LogIO::POST;
1976 0 : for(uInt scale=0; scale < userScaleSizes.nelements(); scale++) {
1977 0 : os << "scale " << scale+1 << " = " << userScaleSizes(scale)
1978 0 : << " pixels" << LogIO::POST;
1979 : }
1980 0 : scaleSizes_p=userScaleSizes;
1981 0 : cleaner_p->setscales(userScaleSizes);
1982 0 : scalesValid_p=true;
1983 :
1984 : } else {
1985 : os << LogIO::SEVERE << "Unknown scale setting algorithm: "
1986 0 : << scaleMethod << LogIO::POST;
1987 0 : return false;
1988 : }
1989 :
1990 0 : return true;
1991 0 : }
1992 :
1993 0 : Bool Deconvolver::clone(const String& imageName, const String& newImageName)
1994 : {
1995 0 : if(!valid()) return false;
1996 0 : LogIO os(LogOrigin("Deconvolver", "clone()", WHERE));
1997 : try {
1998 0 : PagedImage<Float> oldImage(imageName);
1999 0 : PagedImage<Float> newImage(oldImage.shape(), oldImage.coordinates(),
2000 0 : newImageName);
2001 0 : } catch (AipsError x) {
2002 0 : os << LogIO::SEVERE << "Exception: " << x.getMesg() << LogIO::POST;
2003 0 : return false;
2004 0 : }
2005 0 : return true;
2006 0 : }
2007 :
2008 :
2009 :
2010 0 : Bool Deconvolver::convolve(const String& convolvedName,
2011 : const String& modelName)
2012 : {
2013 0 : PagedImage<Float> model(modelName);
2014 0 : PagedImage<Float> convolved(model.shape(),
2015 : model.coordinates(),
2016 0 : convolvedName);
2017 0 : convolver_p->linear( convolved, model );
2018 0 : return true;
2019 0 : };
2020 :
2021 0 : Bool Deconvolver::makegaussian(const String& gaussianName, GaussianBeam& mbeam, Bool normalizeVolume)
2022 : {
2023 0 : LogIO os(LogOrigin("Deconvolver", "makegaussian()", WHERE));
2024 0 : if(!dirty_p) {
2025 : os << LogIO::SEVERE << "Program logic error: Dirty image pointer dirty_p not yet set"
2026 0 : << LogIO::POST;
2027 0 : return false;
2028 : }
2029 0 : PagedImage<Float> gaussian(dirty_p->shape(),
2030 0 : dirty_p->coordinates(),
2031 0 : gaussianName);
2032 0 : gaussian.set(0.0);
2033 0 : gaussian.setUnits(Unit("Jy/pixel"));
2034 0 : putGaussian(gaussian, mbeam);
2035 0 : if(!normalizeVolume){
2036 0 : Float maxpsf=max(gaussian).getFloat();
2037 0 : gaussian.copyData((LatticeExpr<Float>)(gaussian/maxpsf));
2038 : }
2039 : //uInt naxis=gaussian.shape().nelements();
2040 : // StokesImageUnil::Convolve requires an image with four axes
2041 : /*
2042 : if(naxis==2){
2043 : IPosition center = gaussian.shape()/2;
2044 : gaussian.putAt(1.0, center);
2045 : }
2046 : else if(naxis==3){
2047 : IPosition center(3, Int((nx_p/4)*2), Int((ny_p/4)*2),0);
2048 : for (Int k=0; k < gaussian.shape()(2); ++k){
2049 : center(2) = k;
2050 : gaussian.putAt(1.0, center);
2051 : }
2052 : }
2053 : */
2054 : /* if(naxis<4){
2055 : os << LogIO::SEVERE << "Input dirty image: naxis=" <<naxis
2056 : << ". Four axes are required."<< LogIO::POST;
2057 : return false;
2058 :
2059 : }
2060 : else if(naxis==4){
2061 : IPosition center(4, Int((nx_p/4)*2), Int((ny_p/4)*2),0,0);
2062 : for (Int k=0; k < gaussian.shape()(2); ++k){
2063 :
2064 : center(2) = k;
2065 : for(Int j=0; j < gaussian.shape()(3); ++j){
2066 : center(3)=j;
2067 : gaussian.putAt(1.0, center);
2068 : }
2069 : }
2070 :
2071 :
2072 : }
2073 : StokesImageUtil::Convolve(gaussian, mbeam, normalizeVolume);
2074 : */
2075 0 : return true;
2076 0 : };
2077 :
2078 0 : Bool Deconvolver::putGaussian(ImageInterface<Float>& im, const GaussianBeam& beam){
2079 0 : CoordinateSystem cs=im.coordinates();
2080 :
2081 0 : Vector<Int> dirAxes=CoordinateUtil::findDirectionAxes(cs);
2082 0 : DirectionCoordinate dirCoord=cs.directionCoordinate(cs.findCoordinate(Coordinate::DIRECTION));
2083 0 : Vector<Double> cenpix(2, Double(nx_p)/2.0);
2084 0 : cenpix(1)=Double(ny_p)/2.0;
2085 0 : MDirection centre;
2086 0 : dirCoord.toWorld(centre, cenpix);
2087 0 : MVAngle mvRA=centre.getAngle().getValue()(0);
2088 0 : MVAngle mvDEC=centre.getAngle().getValue()(1);
2089 : //cerr << "centre " << cenpix << " " << centre.getRefString() << " " << mvRA(0.0).string(MVAngle::TIME,8) << " " << mvDEC(0.0).string(MVAngle::ANGLE_CLEAN,8) << endl;
2090 : //cerr << "maj min pa " << beam.getMajor() << " " << beam.getMinor() << " " << beam.getPA() << endl;
2091 0 : GaussianShape gshp(centre, beam.getMajor(), beam.getMinor(), beam.getPA());
2092 0 : SkyComponent gcomp(Flux<Double>(1.0), gshp, ConstantSpectrum());
2093 0 : ComponentList cl;
2094 0 : cl.add(gcomp);
2095 0 : ComponentImager::project(im, cl);
2096 0 : return true;
2097 :
2098 :
2099 :
2100 0 : }
2101 :
2102 :
2103 0 : Bool Deconvolver::detached() const
2104 : {
2105 0 : if (dirty_p == 0) {
2106 0 : LogIO os(LogOrigin("Deconvolver", "detached()", WHERE));
2107 : os << LogIO::SEVERE <<
2108 : "Deconvolver is detached - cannot perform operation." << endl <<
2109 0 : "Call Deconvolver.open('dirtyname', 'psfname') to reattach." << LogIO::POST;
2110 0 : return true;
2111 0 : }
2112 0 : return false;
2113 : }
2114 :
2115 0 : Bool Deconvolver::removeTable(const String& tablename) {
2116 :
2117 0 : LogIO os(LogOrigin("Deconvolver", "removeTable()", WHERE));
2118 :
2119 0 : if(Table::isReadable(tablename)) {
2120 0 : if (! Table::isWritable(tablename)) {
2121 : os << LogIO::SEVERE << "Table " << tablename
2122 0 : << " is not writable!: cannot alter it" << LogIO::POST;
2123 0 : return false;
2124 : }
2125 : else {
2126 0 : if (Table::isOpened(tablename)) {
2127 : os << LogIO::SEVERE << "Table " << tablename
2128 : << " is already open in the process. It needs to be closed first"
2129 0 : << LogIO::POST;
2130 0 : return false;
2131 : } else {
2132 0 : Table table(tablename, Table::Update);
2133 0 : if (table.isMultiUsed()) {
2134 : os << LogIO::SEVERE << "Table " << tablename
2135 : << " is already open in another process. It needs to be closed first"
2136 0 : << LogIO::POST;
2137 0 : return false;
2138 : } else {
2139 0 : Table table(tablename, Table::Delete);
2140 0 : }
2141 0 : }
2142 : }
2143 : }
2144 0 : return true;
2145 0 : }
2146 :
2147 0 : Bool Deconvolver::valid() const {
2148 0 : LogIO os(LogOrigin("Deconvolver", "if(!valid()) return false", WHERE));
2149 0 : if(!dirty_p) {
2150 : os << LogIO::SEVERE << "Program logic error: Dirty image pointer dirty_p not yet set"
2151 0 : << LogIO::POST;
2152 0 : return false;
2153 : }
2154 0 : if(!psf_p) {
2155 : os << LogIO::SEVERE << "Program logic error: PSF pointer psf_p not yet set"
2156 0 : << LogIO::POST;
2157 0 : return false;
2158 : }
2159 0 : return true;
2160 0 : }
2161 :
2162 0 : void Deconvolver::findAxes(){
2163 :
2164 0 : CoordinateSystem coordsys= dirty_p->coordinates();
2165 0 : polAxis_p=coordsys.findCoordinate(Coordinate::STOKES);
2166 0 : if(polAxis_p > coordsys.findCoordinate(Coordinate::DIRECTION))
2167 0 : polAxis_p+=1;
2168 0 : chanAxis_p=coordsys.findCoordinate(Coordinate::SPECTRAL);
2169 0 : if(chanAxis_p > coordsys.findCoordinate(Coordinate::DIRECTION))
2170 0 : chanAxis_p+=1;
2171 :
2172 0 : }
2173 :
2174 0 : void Deconvolver::checkMask(ImageInterface<Float>& maskImage, Int& xbeg, Int& xend,
2175 : Int& ybeg, Int& yend){
2176 :
2177 :
2178 0 : LogIO os(LogOrigin("Deconvolver","checkMask",WHERE));
2179 :
2180 0 : xbeg=nx_p/4;
2181 0 : ybeg=ny_p/4;
2182 :
2183 0 : xend=xbeg+nx_p/2-1;
2184 0 : yend=ybeg+ny_p/2-1;
2185 0 : Slicer sl;
2186 0 : if(nchan_p >= 1){
2187 :
2188 0 : if(npol_p > 0 ){
2189 0 : IPosition blc(4,0,0,0,0);
2190 0 : blc(chanAxis_p)=0;
2191 0 : blc(polAxis_p)=0;
2192 0 : IPosition trc(4, nx_p-1, ny_p-1, 0, 0);
2193 0 : trc(chanAxis_p)=0;
2194 0 : trc(polAxis_p)=0;
2195 0 : sl=Slicer(blc, trc, Slicer::endIsLast);
2196 0 : }
2197 : else{
2198 0 : IPosition blc(3, 0, 0, 0);
2199 0 : IPosition trc(3, nx_p-1, ny_p-1, 0);
2200 0 : sl=Slicer(blc, trc, Slicer::endIsLast);
2201 0 : }
2202 : }
2203 : else{
2204 0 : if(npol_p > 0 ){
2205 0 : IPosition blc(3, 0, 0, 0);
2206 0 : IPosition trc(3, nx_p-1, ny_p-1, 0);
2207 0 : sl=Slicer(blc, trc, Slicer::endIsLast);
2208 0 : }
2209 : else{
2210 0 : IPosition blc(2, 0, 0);
2211 0 : IPosition trc(2, nx_p-1, ny_p-1);
2212 0 : }
2213 : }
2214 :
2215 :
2216 0 : Matrix<Float> mask= maskImage.getSlice(sl, true);
2217 : // ignore mask if none exists
2218 0 : if(max(mask) < 0.000001) {
2219 : os << "Mask seems to be empty; will CLEAN inner quarter"
2220 0 : << LogIO::WARN;
2221 0 : return;
2222 : }
2223 : // Now read the mask and determine the bounding box
2224 :
2225 0 : xbeg=nx_p-1;
2226 0 : ybeg=ny_p-1;
2227 0 : xend=0;
2228 0 : yend=0;
2229 :
2230 :
2231 0 : for (Int iy=0;iy<ny_p;iy++) {
2232 0 : for (Int ix=0;ix<nx_p;ix++) {
2233 0 : if(mask(ix,iy)>0.000001) {
2234 0 : xbeg=min(xbeg,ix);
2235 0 : ybeg=min(ybeg,iy);
2236 0 : xend=max(xend,ix);
2237 0 : yend=max(yend,iy);
2238 :
2239 : }
2240 : }
2241 : }
2242 : // Now have possible BLC. Make sure that we don't go over the
2243 : // edge later
2244 0 : if(((xend - xbeg)>nx_p/2) && (!fullPlane_p)) {
2245 0 : xbeg=nx_p/4-1; //if larger than quarter take inner of mask
2246 0 : os << LogIO::WARN << "Mask span over more than half the x-axis: Considering inner half of the x-axis" << LogIO::POST;
2247 : }
2248 0 : if(((yend - ybeg)>ny_p/2) && (!fullPlane_p)) {
2249 0 : ybeg=ny_p/4-1;
2250 0 : os << LogIO::WARN << "Mask span over more than half the y-axis: Considering inner half of the y-axis" << LogIO::POST;
2251 : }
2252 :
2253 : // Just making sure we are within limits...
2254 0 : if(fullPlane_p){
2255 0 : xend=min(xend,nx_p-1);
2256 0 : yend=min(yend,ny_p-1);
2257 : }
2258 : else{
2259 0 : xend=min(xend,xbeg+nx_p/2-1);
2260 0 : yend=min(yend,ybeg+ny_p/2-1);
2261 : }
2262 :
2263 :
2264 :
2265 0 : }
2266 :
2267 : ///////////////////////////////////////////////////////////////////////////////////////////////////////////////
2268 : // Multi-Term Clean algorithm with Taylor-Polynomial basis functions
2269 : ///////////////////////////////////////////////////////////////////////////////////////////////////////////////
2270 0 : Bool Deconvolver::mtopen(const Int nTaylor,
2271 : const Vector<Float>& userScaleSizes,
2272 : const Vector<String>& psfs)
2273 : {
2274 :
2275 : // if(!valid()) return false; //make MT version
2276 0 : LogIO os(LogOrigin("Deconvolver", "mtopen()", WHERE));
2277 :
2278 : // Check for already-open mt-deconvolver
2279 :
2280 0 : if(mt_nterms_p != -1)
2281 : {
2282 0 : os << LogIO::WARN << "Multi-term Deconvolver is already open, and set-up for " << mt_nterms_p << " Taylor-terms. Please close the deconvolver and reopen, or continue with cleaning." << LogIO::POST;
2283 0 : return false;
2284 : }
2285 :
2286 : // Check for valid ntaylor
2287 0 : if( nTaylor <=0 )
2288 : {
2289 0 : os << LogIO::SEVERE << "nTaylor must be at-least 1" << LogIO::POST;
2290 0 : return false;
2291 : }
2292 :
2293 0 : mt_nterms_p = nTaylor;
2294 0 : uInt mt_npsftaylor = 2*nTaylor-1;
2295 :
2296 : // Check if the correct number of PSFs exist.
2297 0 : if( psfs.nelements() != mt_npsftaylor )
2298 : {
2299 0 : os << LogIO::SEVERE << "For " << mt_nterms_p << " Taylor terms, " << mt_npsftaylor << " PSFs are needed to populate the Hessian matrix. " << LogIO::POST;
2300 0 : return false;
2301 : }
2302 :
2303 0 : os << "Initializing MT-Clean with " << mt_nterms_p << " Taylor terms, " << userScaleSizes.nelements() << " scales and " << psfs.nelements() << " unique Hessian elements (psfs)" << LogIO::POST;
2304 :
2305 : // Open all the PSFs...
2306 0 : Block<CountedPtr<PagedImage<Float> > > mt_psfs(mt_npsftaylor);
2307 :
2308 : // Open the first PSF and extract shape info.
2309 0 : mt_psfs[0] = new PagedImage<Float>(psfs[0]);
2310 0 : AlwaysAssert(&*mt_psfs[0], AipsError);
2311 0 : nx_p=mt_psfs[0]->shape()(0);
2312 0 : ny_p=mt_psfs[0]->shape()(1);
2313 : //mt_psfs_p[0].setMaximumCacheSize(2*nx_p*ny_p);
2314 : //mt_psfs_p[0].setCacheSizeInTiles(10000);
2315 :
2316 : // Open the other PSFs and verify shape consistency
2317 0 : for(uInt i=1; i<mt_npsftaylor;i++)
2318 : {
2319 0 : mt_psfs[i] = new PagedImage<Float>(psfs[i]);
2320 0 : AlwaysAssert(&*mt_psfs[i], AipsError);
2321 0 : if( mt_psfs[i]->shape()(0) != nx_p || mt_psfs[i]->shape()(1) != ny_p )
2322 : {
2323 0 : os << LogIO::SEVERE << "Supplied PSFs are of different shapes. Please try again." << LogIO::POST;
2324 0 : mt_psfs.resize(0); // Not sure if this is safe, with CountedPtrs. WARN.
2325 0 : return false;
2326 : }
2327 : }
2328 :
2329 : // Initialize the Multi-Term Cleaner
2330 : try
2331 : {
2332 0 : mt_cleaner_p.setntaylorterms(mt_nterms_p);
2333 0 : mt_cleaner_p.setscales(userScaleSizes);
2334 0 : mt_cleaner_p.initialise(nx_p,ny_p); // allocates memory once....
2335 : }
2336 0 : catch (AipsError x)
2337 : {
2338 : os << LogIO::WARN << "Cannot allocate required memory for Multi-Term minor cycle"
2339 0 : << LogIO::POST;
2340 0 : return false;
2341 0 : }
2342 :
2343 : // Send the PSFs into the Multi-Term Matrix Cleaner
2344 0 : for (uInt order=0;order<mt_npsftaylor;order++)
2345 : {
2346 0 : Matrix<Float> tempMat;
2347 0 : Array<Float> tempArr;
2348 0 : (mt_psfs[order])->get(tempArr,true); // by reference.
2349 0 : tempMat.reference(tempArr);
2350 :
2351 0 : mt_cleaner_p.setpsf( order , tempMat );
2352 0 : }
2353 :
2354 : // Compute Hessian elements for Taylor terms and Scales ( convolutions ), take the Hessian block-diagonal approximation and inverse, Check for invertability.
2355 0 : Int ret = mt_cleaner_p.computeHessianPeak();
2356 0 : if (ret != 0)
2357 : {
2358 0 : os << LogIO::SEVERE << "Cannot Invert Hessian matrix. Please close the deconvolver and supply different PSFs." << LogIO::POST;
2359 0 : return false;
2360 : }
2361 :
2362 : // mt_psfs goes out of scope and CountedPtrs delete themselves automatically.
2363 :
2364 0 : mt_valid_p=true;
2365 :
2366 0 : return true;
2367 0 : }
2368 :
2369 0 : Bool Deconvolver::mtclean(const Vector<String>& residuals,
2370 : const Vector<String>& models,
2371 : const Int niter,
2372 : const Float gain,
2373 : const Quantity& threshold,
2374 : const Bool /*displayProgress*/,
2375 : const String& mask,
2376 : Float& maxResidual, Int& iterationsDone)
2377 : {
2378 :
2379 0 : LogIO os(LogOrigin("Deconvolver", "mtclean()", WHERE));
2380 :
2381 0 : if(mt_valid_p==false)
2382 : {
2383 0 : os << LogIO::WARN << "Multi-Term Deconvolver is not initialized yet. Please call 'mtopen()' first" << LogIO::POST;
2384 0 : return false;
2385 : }
2386 :
2387 0 : os << "Running MT-Clean with " << mt_nterms_p << " Taylor terms " << LogIO::POST;
2388 :
2389 : // Send in the mask.
2390 0 : if( mask != String("") )
2391 : {
2392 0 : PagedImage<Float> mt_mask(mask);
2393 0 : Matrix<Float> tempMat;
2394 0 : Array<Float> tempArr;
2395 0 : (mt_mask).get(tempArr,true);
2396 0 : tempMat.reference(tempArr);
2397 :
2398 0 : mt_cleaner_p.setmask( tempMat );
2399 0 : }
2400 :
2401 : // Send in the residuals and model images.
2402 0 : Block<CountedPtr<PagedImage<Float> > > mt_residuals(mt_nterms_p);
2403 0 : Block<CountedPtr<PagedImage<Float> > > mt_models(mt_nterms_p);
2404 0 : for (Int i=0;i<mt_nterms_p;i++)
2405 : {
2406 0 : mt_residuals[i] = new PagedImage<Float>(residuals[i]);
2407 0 : AlwaysAssert(&*mt_residuals[i], AipsError);
2408 0 : if( mt_residuals[i]->shape()(0) != nx_p || mt_residuals[i]->shape()(1) != ny_p )
2409 : {
2410 0 : os << LogIO::SEVERE << "Supplied Residual images don't match PSF shapes." << LogIO::POST;
2411 0 : mt_residuals.resize(0); // Not sure if this is safe, with CountedPtrs. WARN.
2412 0 : return false;
2413 : }
2414 :
2415 0 : mt_models[i] = new PagedImage<Float>(models[i]);
2416 0 : AlwaysAssert(&*mt_models[i], AipsError);
2417 0 : if( mt_models[i]->shape()(0) != nx_p || mt_models[i]->shape()(1) != ny_p )
2418 : {
2419 0 : os << LogIO::SEVERE << "Supplied Model images don't match PSF shapes." << LogIO::POST;
2420 0 : mt_models.resize(0); // Not sure if this is safe, with CountedPtrs. WARN.
2421 0 : return false;
2422 : }
2423 :
2424 :
2425 0 : Matrix<Float> tempMat;
2426 0 : Array<Float> tempArr;
2427 :
2428 0 : (mt_residuals[i])->get(tempArr,true); // by reference.
2429 0 : tempMat.reference(tempArr);
2430 0 : mt_cleaner_p.setresidual( i , tempMat );
2431 :
2432 0 : (mt_models[i])->get(tempArr,true); // by reference.
2433 0 : tempMat.reference(tempArr);
2434 0 : mt_cleaner_p.setmodel( i , tempMat );
2435 :
2436 0 : } // end of for i
2437 :
2438 :
2439 : // Fill in maxResidual and iterationsDone
2440 :
2441 0 : Float fractionOfPsf=0.05;
2442 :
2443 : // Call the cleaner
2444 0 : iterationsDone = mt_cleaner_p.mtclean(niter, fractionOfPsf, gain, threshold.getValue(String("Jy")));
2445 :
2446 : // Get back the model (this is not held in MTMC by reference, because of 'incremental=T/F' logic....
2447 0 : for (Int order=0;order<mt_nterms_p;order++)
2448 : {
2449 0 : Matrix<Float> tempMod;
2450 0 : mt_cleaner_p.getmodel(order,tempMod);
2451 0 : mt_models[order]->put(tempMod);
2452 :
2453 : // Also get the updated residuals. Not really needed, but good to have.
2454 0 : mt_cleaner_p.getresidual(order,tempMod);
2455 0 : mt_residuals[order]->put(tempMod);
2456 :
2457 0 : if(order==0) maxResidual = max(fabs(tempMod));
2458 :
2459 0 : }
2460 :
2461 0 : return maxResidual <= threshold.getValue(String("Jy"));
2462 0 : }
2463 :
2464 :
2465 0 : Bool Deconvolver::mtrestore(const Vector<String>& models,
2466 : const Vector<String>& residuals,
2467 : const Vector<String>& images,
2468 : GaussianBeam& mbeam)
2469 : {
2470 :
2471 0 : LogIO os(LogOrigin("Deconvolver", "mtrestore()", WHERE));
2472 :
2473 : // check that the mt_cleaner_p is alive..
2474 0 : if(mt_valid_p==false)
2475 : {
2476 0 : os << LogIO::WARN << "Multi-Term Deconvolver is not initialized yet. Please call 'mtopen()' first" << LogIO::POST;
2477 0 : return false;
2478 : }
2479 :
2480 0 : os << "Restoring Taylor-coefficient images" << LogIO::POST;
2481 :
2482 0 : Block<CountedPtr<PagedImage<Float> > > mt_residuals(mt_nterms_p);
2483 0 : for (Int i=0;i<mt_nterms_p;i++)
2484 : {
2485 0 : mt_residuals[i] = new PagedImage<Float>(residuals[i]);
2486 0 : AlwaysAssert(&*mt_residuals[i], AipsError);
2487 0 : if( mt_residuals[i]->shape()(0) != nx_p || mt_residuals[i]->shape()(1) != ny_p )
2488 : {
2489 0 : os << LogIO::SEVERE << "Supplied Residual images don't match PSF shapes." << LogIO::POST;
2490 0 : mt_residuals.resize(0); // Not sure if this is safe, with CountedPtrs. WARN.
2491 0 : return false;
2492 : }
2493 : }
2494 :
2495 : // Calculate principal solution on the residuals.
2496 : //// Send in new residuals
2497 0 : for (Int order=0;order<mt_nterms_p;order++)
2498 : {
2499 0 : Matrix<Float> tempMat;
2500 0 : Array<Float> tempArr;
2501 0 : (mt_residuals[order])->get(tempArr,true); // by reference.
2502 0 : tempMat.reference(tempArr);
2503 0 : mt_cleaner_p.setresidual( order , tempMat );
2504 0 : }
2505 :
2506 : //// Compute p-soln
2507 0 : mt_cleaner_p.computeprincipalsolution();
2508 :
2509 : //// Get residuals out
2510 0 : for (Int order=0;order<mt_nterms_p;order++)
2511 : {
2512 0 : Matrix<Float> tempMat;
2513 0 : Matrix<Float> tempMod;
2514 0 : mt_cleaner_p.getresidual(order,tempMod);
2515 0 : mt_residuals[order]->put(tempMod);
2516 0 : }
2517 :
2518 : // Convolve models with the clean beam and add new residuals.... per term
2519 0 : for (Int order=0;order<mt_nterms_p;order++)
2520 : {
2521 0 : PagedImage<Float> thismodel(models[order]);
2522 0 : PagedImage<Float> thisimage( thismodel.shape(), thismodel.coordinates(), images[order]);
2523 :
2524 0 : LatticeExpr<Float> cop(thismodel);
2525 0 : thisimage.copyData(cop);
2526 :
2527 0 : StokesImageUtil::Convolve(thisimage, mbeam);
2528 :
2529 0 : LatticeExpr<Float> le(thisimage+( *mt_residuals[order] ));
2530 0 : thisimage.copyData(le);
2531 :
2532 0 : ImageInfo ii = thisimage.imageInfo();
2533 0 : ii.setRestoringBeam(mbeam);
2534 0 : thisimage.setImageInfo(ii);
2535 0 : thisimage.setUnits(Unit("Jy/beam"));
2536 0 : thisimage.table().unmarkForDelete();
2537 :
2538 0 : }
2539 :
2540 0 : return true;
2541 0 : }
2542 :
2543 : // This code duplicates WBCleanImageSkyModel::calculateAlphaBeta
2544 : // Eventually, This Deconvolver code will replace what's in WBCleanImageSkyModel.cc
2545 : // This function must be callable stand-alone, and must not use private vars.
2546 : // This is to allow easy recalculation of alpha with different thresholds
2547 0 : Bool Deconvolver::mtcalcpowerlaw(const Vector<String>& images,
2548 : const Vector<String>& residuals,
2549 : const String& alphaname,
2550 : const String& betaname,
2551 : const Quantity& threshold,
2552 : const Bool calcerror)
2553 : {
2554 0 : LogIO os(LogOrigin("Deconvolver", "mtcalcpowerlaw()", WHERE));
2555 :
2556 0 : uInt ntaylor = images.nelements();
2557 :
2558 0 : if(ntaylor<2)
2559 : {
2560 0 : os << LogIO::SEVERE << "Please enter at-least two Taylor-coefficient images" << LogIO::POST;
2561 0 : return false;
2562 : }
2563 : else
2564 : {
2565 :
2566 : // Check that the images exist on disk
2567 :
2568 0 : os << "Calculating spectral index";
2569 0 : if(ntaylor>3) os << " and spectral curvature";
2570 0 : os << " from " << ntaylor << " restored images, using a mask defined by a threshold of " << threshold.getValue(String("Jy")) << " Jy " << LogIO::POST;
2571 : }
2572 :
2573 : // Open restored images
2574 0 : PagedImage<Float> imtaylor0(images[0]);
2575 0 : PagedImage<Float> imtaylor1(images[1]);
2576 :
2577 : // Check shapes
2578 0 : if( imtaylor0.shape() != imtaylor1.shape() )
2579 : {
2580 0 : os << LogIO::SEVERE << "Taylor-coefficient image shapes must match." << LogIO::POST;
2581 0 : return false;
2582 : }
2583 :
2584 : // Create empty alpha image
2585 0 : PagedImage<Float> imalpha(imtaylor0.shape(),imtaylor0.coordinates(),alphaname);
2586 0 : imalpha.set(0.0);
2587 :
2588 0 : Float specthreshold = threshold.getValue(String("Jy"));
2589 :
2590 : // Create a mask - make this adapt to the signal-to-noise
2591 0 : LatticeExpr<Float> mask1(iif((imtaylor0)>(specthreshold),1.0,0.0));
2592 0 : LatticeExpr<Float> mask0(iif((imtaylor0)>(specthreshold),0.0,1.0));
2593 :
2594 : /////// Calculate alpha
2595 0 : LatticeExpr<Float> alphacalc( ((imtaylor1)*mask1)/((imtaylor0)+(mask0)) );
2596 0 : imalpha.copyData(alphacalc);
2597 :
2598 : // Set the restoring beam for alpha
2599 0 : ImageInfo ii = imalpha.imageInfo();
2600 0 : ii.setRestoringBeam( (imtaylor0.imageInfo()).restoringBeam() );
2601 0 : imalpha.setImageInfo(ii);
2602 : //imalpha.setUnits(Unit("Spectral Index"));
2603 0 : imalpha.table().unmarkForDelete();
2604 :
2605 : // Make a mask for the alpha image
2606 0 : LatticeExpr<Bool> lemask(iif((imtaylor0 > specthreshold) , true, false));
2607 :
2608 0 : createMask(lemask, imalpha);
2609 0 : os << "Written Spectral Index Image : " << alphaname << LogIO::POST;
2610 :
2611 :
2612 : ////// Calculate error on alpha, if requested.
2613 0 : if(calcerror)
2614 : {
2615 0 : if( residuals.nelements() != ntaylor )
2616 : {
2617 0 : os << LogIO::WARN << "Number of residual images is different from number of restored images. Not calculating alpha-error map." << LogIO::POST;
2618 : }
2619 : else
2620 : {
2621 0 : PagedImage<Float> imalphaerror(imtaylor0.shape(),imtaylor0.coordinates(),alphaname+String(".error"));
2622 0 : imalphaerror.set(0.0);
2623 :
2624 : /* Open residual images */
2625 0 : PagedImage<Float> residual0(residuals[0]);
2626 0 : PagedImage<Float> residual1(residuals[1]);
2627 :
2628 0 : if( residual0.shape() != residual1.shape() || residual0.shape() != imtaylor0.shape() )
2629 : {
2630 0 : os << LogIO::WARN << "Shapes of residual images (and/or restored images) do not match. Not calculating alpha-error map." << LogIO::POST;
2631 : }
2632 : else
2633 : {
2634 0 : LatticeExpr<Float> alphacalcerror( abs(alphacalc) * sqrt( ( (residual0*mask1)/(imtaylor0+mask0) )*( (residual0*mask1)/(imtaylor0+mask0) ) + ( (residual1*mask1)/(imtaylor1+mask0) )*( (residual1*mask1)/(imtaylor1+mask0) ) ) );
2635 0 : imalphaerror.copyData(alphacalcerror);
2636 0 : imalphaerror.setImageInfo(ii);
2637 0 : createMask(lemask, imalphaerror);
2638 0 : imalphaerror.table().unmarkForDelete();
2639 0 : os << "Written Spectral Index Error Image : " << alphaname << ".error" << LogIO::POST;
2640 :
2641 : // mergeDataError( imalpha, imalphaerror, alphaerrorname+".new" );
2642 0 : }
2643 :
2644 0 : }
2645 :
2646 : }// end if calcerror
2647 :
2648 : ////// Calculate beta, if enough images are given.
2649 0 : if(ntaylor>2)
2650 : {
2651 0 : PagedImage<Float> imbeta(imtaylor0.shape(),imtaylor0.coordinates(),betaname);
2652 0 : imbeta.set(0.0);
2653 0 : PagedImage<Float> imtaylor2(images[2]);
2654 0 : if( imtaylor2.shape() != imtaylor0.shape() )
2655 : {
2656 0 : os << LogIO::WARN << "Restored image shapes do not match. Not calculating 'beta'" << LogIO::POST;
2657 : }
2658 : else
2659 : {
2660 0 : LatticeExpr<Float> betacalc( ((imtaylor2)*mask1)/((imtaylor0)+(mask0))-0.5*(imalpha)*(imalpha-1.0) );
2661 0 : imbeta.copyData(betacalc);
2662 0 : imbeta.setImageInfo(ii);
2663 : //imbeta.setUnits(Unit("Spectral Curvature"));
2664 0 : createMask(lemask, imbeta);
2665 0 : imbeta.table().unmarkForDelete();
2666 :
2667 0 : os << "Written Spectral Curvature Image : " << betaname << LogIO::POST;
2668 0 : }
2669 0 : }
2670 :
2671 0 : return true;
2672 0 : }
2673 :
2674 : // This is also a copy of WBCleanImageSkyModel::createMask
2675 : // Eventually, WBCleanImageSkyModel must use this Deconvolver version.
2676 0 : Bool Deconvolver::createMask(LatticeExpr<Bool> &lemask, ImageInterface<Float> &outimage)
2677 : {
2678 0 : ImageRegion outreg = outimage.makeMask("mask0",false,true);
2679 0 : LCRegion& outmask=outreg.asMask();
2680 0 : outmask.copyData(lemask);
2681 0 : outimage.defineRegion("mask0",outreg, RegionHandler::Masks, true);
2682 0 : outimage.setDefaultMask("mask0");
2683 0 : return true;
2684 0 : }
2685 :
2686 :
2687 : } //# NAMESPACE CASA - END
|
