Line data Source code
1 : //# CEMemImageSkyModel.cc: Implementation of CEMemImageSkyModel class
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: aips2-request@nrao.edu.
21 : //# Postal address: AIPS++ Project Office
22 : //# National Radio Astronomy Observatory
23 : //# 520 Edgemont Road
24 : //# Charlottesville, VA 22903-2475 USA
25 : //#
26 : //# $Id$
27 :
28 : #include <casacore/casa/Arrays/ArrayMath.h>
29 : #include <synthesis/MeasurementComponents/CEMemImageSkyModel.h>
30 : #include <casacore/images/Images/PagedImage.h>
31 : #include <casacore/casa/OS/File.h>
32 : #include <casacore/lattices/Lattices/LatticeStepper.h>
33 : #include <casacore/lattices/Lattices/LatticeIterator.h>
34 : #include <casacore/lattices/LEL/LatticeExpr.h>
35 : #include <casacore/lattices/LEL/LatticeExprNode.h>
36 : #include <casacore/lattices/LRegions/LCBox.h>
37 : #include <casacore/lattices/Lattices/SubLattice.h>
38 : #include <synthesis/MeasurementEquations/SkyEquation.h>
39 : #include <casacore/casa/Exceptions/Error.h>
40 : #include <casacore/casa/BasicSL/String.h>
41 : #include <casacore/casa/Utilities/Assert.h>
42 :
43 : #include <sstream>
44 :
45 : #include <casacore/casa/Logging/LogMessage.h>
46 : #include <casacore/casa/Logging/LogIO.h>
47 : #include <casacore/casa/Logging/LogSink.h>
48 :
49 : #include <synthesis/MeasurementEquations/LatConvEquation.h>
50 : #include <synthesis/MeasurementEquations/IncCEMemModel.h>
51 : #include <synthesis/MeasurementEquations/CEMemProgress.h>
52 :
53 :
54 : using namespace casacore;
55 : namespace casa { //# NAMESPACE CASA - BEGIN
56 :
57 0 : CEMemImageSkyModel::
58 : CEMemImageSkyModel(Float sigma,
59 : Float targetFlux,
60 : Bool constrainFlux,
61 : const Vector<String>& priors,
62 0 : const String& entropy)
63 : :
64 : itsSigma(sigma),
65 : itsTargetFlux(targetFlux),
66 : itsConstrainFlux(constrainFlux),
67 : itsPrior(priors),
68 : itsEntropy(entropy),
69 : itsInitializeModel(true),
70 0 : itsProgress(0)
71 : {
72 0 : };
73 :
74 0 : CEMemImageSkyModel::~CEMemImageSkyModel()
75 : {
76 0 : }
77 :
78 : // Mem solver
79 0 : Bool CEMemImageSkyModel::solve(SkyEquation& se) {
80 :
81 0 : LogIO os(LogOrigin("CEMemImageSkyModel","solve",WHERE));
82 :
83 :
84 0 : if(numberOfModels()>1) {
85 0 : os << "Cannot process more than one field" << LogIO::EXCEPTION;
86 : }
87 :
88 : // Make the residual image
89 0 : makeNewtonRaphsonStep(se);
90 :
91 : //Make the PSF
92 0 : makeApproxPSFs(se);
93 :
94 0 : Int nx=image(0).shape()(0);
95 0 : Int ny=image(0).shape()(1);
96 0 : Int npol=image(0).shape()(2);
97 0 : Int nchan=image(0).shape()(3);
98 :
99 : PagedImage<Float>* priorImagePtr;
100 0 : priorImagePtr=0;
101 0 : if(itsPrior.nelements()>0) {
102 0 : if(Table::isReadable(itsPrior(0)))
103 0 : priorImagePtr=new PagedImage<Float>(itsPrior(0));
104 : }
105 :
106 0 : AlwaysAssert((npol==1)||(npol==2)||(npol==4), AipsError);
107 :
108 : // Now read the mask and determine the bounding box
109 0 : Int xbeg=nx/4;
110 0 : Int ybeg=ny/4;
111 :
112 0 : Int xend=xbeg+nx/2-1;
113 0 : Int yend=ybeg+ny/2-1;
114 :
115 0 : if(hasMask(0)) {
116 0 : AlwaysAssert(mask(0).shape()(0)==nx, AipsError);
117 0 : AlwaysAssert(mask(0).shape()(1)==ny, AipsError);
118 :
119 0 : LatticeStepper mls(mask(0).shape(),
120 0 : IPosition(4, nx, ny, 1, 1),
121 0 : IPosition(4, 0, 1, 3, 2));
122 0 : RO_LatticeIterator<Float> maskli(mask(0), mls);
123 0 : maskli.reset();
124 0 : xbeg=nx-1;
125 0 : ybeg=ny-1;
126 0 : xend=0;
127 0 : yend=0;
128 0 : for (Int iy=0;iy<ny;iy++) {
129 0 : for (Int ix=0;ix<nx;ix++) {
130 0 : if(maskli.matrixCursor()(ix,iy)>0.000001) {
131 0 : xbeg=min(xbeg,ix);
132 0 : ybeg=min(ybeg,iy);
133 0 : xend=max(xend,ix);
134 0 : yend=max(yend,iy);
135 : }
136 : }
137 : }
138 : // Now have possible BLC. Make sure that we don't go over the
139 : // edge later
140 0 : if((xend - xbeg)>nx/2) {
141 0 : xbeg=nx/4-1; //if larger than quarter take inner of mask
142 0 : os << LogIO::WARN << "Mask span over more than half the x-axis: Considering inner half of the x-axis" << LogIO::POST;
143 : }
144 0 : if((yend - ybeg)>ny/2) {
145 0 : ybeg=ny/4-1;
146 0 : os << LogIO::WARN << "Mask span over more than half the y-axis: Considering inner half of the y-axis" << LogIO::POST;
147 : }
148 0 : xend=min(xend,xbeg+nx/2-1);
149 0 : yend=min(yend,ybeg+ny/2-1);
150 : }
151 :
152 : // Mask logic used here: one mask for all Stokes, for all Channels
153 0 : SubLattice<Float>* mask_sl_p = 0;
154 0 : if (hasMask(0) & (xend > xbeg) && (yend > ybeg) ) {
155 0 : LCBox maskbox (IPosition(4, xbeg, ybeg, 0, 0),
156 0 : IPosition(4, xend, yend, 0, 0), mask(0).shape());
157 0 : mask_sl_p = new SubLattice<Float> (mask(0), maskbox, false);
158 : }
159 :
160 :
161 0 : for (Int ichan=0; ichan < nchan; ichan++) {
162 0 : LCBox imagebox(IPosition(4, xbeg, ybeg, 0, ichan),
163 0 : IPosition(4, xend, yend, npol-1, ichan),
164 0 : image(0).shape());
165 0 : LCBox psfbox(IPosition(4, 0, 0, 0, ichan),
166 0 : IPosition(4, nx-1, ny-1, 0, ichan),
167 0 : PSF(0).shape());
168 :
169 0 : SubLattice<Float> psf_sl (PSF(0), psfbox, false);
170 0 : SubLattice<Float> residual_sl (residual(0), imagebox, true);
171 0 : SubLattice<Float> model_sl (image(0), imagebox, true);
172 :
173 0 : TempLattice<Float> dirty_sl( residual_sl.shape());
174 0 : dirty_sl.copyData(residual_sl);
175 :
176 : ResidualEquation<Lattice<Float> > * eqn_p;
177 : Float psfmax;
178 : {
179 0 : LatticeExprNode node = max(psf_sl);
180 0 : psfmax = node.getFloat();
181 : }
182 0 : if(nchan>1) {
183 0 : os<<"Processing channel "<<ichan+1<<" of "<<nchan<<LogIO::POST;
184 : }
185 0 : if(psfmax==0.0) {
186 0 : os << "No data for this channel: skipping" << LogIO::POST;
187 : } else {
188 0 : eqn_p = new LatConvEquation (psf_sl, dirty_sl);
189 :
190 : IncEntropy * myEnt_p;
191 0 : String entString = entropy();
192 0 : if(entString=="entropy") {
193 : os << "Deconvolving image using maximum entropy algorithm"
194 0 : << LogIO::POST;
195 0 : myEnt_p = new IncEntropyI;
196 : }
197 0 : else if (entString=="emptiness") {
198 0 : myEnt_p = new IncEntropyEmptiness;
199 : }
200 : else {
201 0 : os << " Known MEM entropies: entropy | emptiness " << LogIO::POST;
202 : os << LogIO::SEVERE << "Unknown MEM entropy: " << entString
203 0 : << LogIO::POST;
204 0 : return false;
205 : }
206 :
207 0 : TempLattice<Float> zero (model_sl.shape());
208 0 : zero.set(0.0);
209 :
210 0 : IncCEMemModel memer(*myEnt_p, zero, model_sl, numberIterations(),
211 : sigma(),
212 0 : targetFlux(), constrainFlux(),
213 0 : initializeModel(), false );
214 :
215 0 : if(priorImagePtr!=0) {
216 0 : memer.setPrior(*priorImagePtr);
217 : }
218 0 : if (mask_sl_p != 0 ) {
219 0 : memer.setMask( *mask_sl_p );
220 : }
221 :
222 0 : if (displayProgress_p) {
223 0 : itsProgress = new CEMemProgress( pgplotter_p );
224 0 : memer.setProgress(*itsProgress);
225 : }
226 :
227 0 : memer.solve(*eqn_p);
228 : // memer.setChoose(false); // not yet implemented!
229 : // os << "Mem used " << cleaner.numberIterations() << " iterations"
230 : // << " to get to a max residual of " << cleaner.threshold()
231 : // << LogIO::POST;
232 :
233 0 : eqn_p->residual(residual_sl, memer);
234 0 : if (itsProgress) {
235 0 : delete itsProgress;
236 : }
237 : }
238 : }
239 0 : if (priorImagePtr) delete priorImagePtr; priorImagePtr=0;
240 0 : if (mask_sl_p) delete mask_sl_p; mask_sl_p=0;
241 0 : return(true);
242 : };
243 :
244 :
245 : } //# NAMESPACE CASA - END
246 :
|
