Line data Source code
1 : //# WFCleanImageSkyModel.cc: Implementation of WFCleanImageSkyModel 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: 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 : //# $Id$
27 :
28 : #include <casacore/casa/Arrays/ArrayMath.h>
29 : #include <synthesis/MeasurementComponents/WFCleanImageSkyModel.h>
30 : //#include <synthesis/MeasurementComponents/MFCleanImageSkyModel.h>
31 : #include <casacore/images/Images/PagedImage.h>
32 : #include <casacore/images/Images/ImageInterface.h>
33 : #include <casacore/casa/OS/File.h>
34 : #include <casacore/images/Images/SubImage.h>
35 : #include <casacore/lattices/Lattices/LatticeStepper.h>
36 : #include <casacore/lattices/Lattices/LatticeIterator.h>
37 : #include <casacore/lattices/LEL/LatticeExpr.h>
38 : #include <casacore/lattices/LRegions/LCBox.h>
39 : #include <casacore/casa/Exceptions/Error.h>
40 : #include <casacore/casa/BasicSL/String.h>
41 : #include <casacore/casa/Utilities/Assert.h>
42 : #include <casacore/coordinates/Coordinates/DirectionCoordinate.h>
43 : #include <casacore/coordinates/Coordinates/CoordinateSystem.h>
44 : #include <synthesis/TransformMachines/StokesImageUtil.h>
45 : #include <sstream>
46 :
47 : #include <casacore/casa/Logging/LogMessage.h>
48 : #include <casacore/casa/Logging/LogSink.h>
49 : #include <casacore/casa/Logging/LogIO.h>
50 :
51 :
52 : using namespace casacore;
53 : namespace casa { //# NAMESPACE CASA - BEGIN
54 :
55 0 : WFCleanImageSkyModel::WFCleanImageSkyModel():
56 0 : MFCleanImageSkyModel(), nfacets_p(1), facets_p(1) /*,
57 : largeMem_p(false)*/ {
58 :
59 0 : imageImage_p=0;
60 0 : residualImage_p=0;
61 0 : maskImage_p=0;
62 0 : };
63 :
64 0 : WFCleanImageSkyModel::WFCleanImageSkyModel(const Int facets, Bool ):
65 0 : MFCleanImageSkyModel(), nfacets_p(facets*facets), facets_p(facets) /*,
66 : largeMem_p(largeMemory) */
67 : {
68 0 : imageImage_p=0;
69 0 : residualImage_p=0;
70 0 : maskImage_p=0;
71 :
72 0 : };
73 :
74 0 : Int WFCleanImageSkyModel::add(ImageInterface<Float>& image,
75 : const Int maxNumXfr)
76 : {
77 : // Is this the first image?
78 0 : if(imageImage_p.null()) {
79 : // Create the facet images and copy the relevant region from the
80 : // original image
81 0 : imageImage_p=CountedPtr<ImageInterface<Float> >(&image, false);
82 0 : facetImages_p.resize(nfacets_p);
83 0 : for (Int facet=0;facet<nfacets_p;facet++) {
84 0 : facetImages_p[facet] = makeFacet(facet, image);
85 0 : AlwaysAssert(!facetImages_p[facet].null(), AipsError);
86 0 : AlwaysAssert(MFCleanImageSkyModel::add(*facetImages_p[facet], maxNumXfr)==facet,
87 : AipsError);
88 : }
89 0 : return 0;
90 : }
91 : else {
92 : // All other images are added as is
93 0 : Int result = MFCleanImageSkyModel::add(image, maxNumXfr);
94 0 : return result-nfacets_p+1;
95 : }
96 : };
97 :
98 0 : Bool WFCleanImageSkyModel::addMask(Int image, ImageInterface<Float>& mask)
99 : {
100 : // Create the facet images and copy the relevant region from the
101 : // original image
102 0 : AlwaysAssert(image>-1, AipsError);
103 0 : if(maskImage_p.null()) {
104 0 : maskImage_p=CountedPtr<ImageInterface<Float> >(&mask, false);
105 0 : facetMaskImages_p.resize(nfacets_p);
106 0 : for (Int facet=0;facet<nfacets_p;facet++) {
107 0 : facetMaskImages_p[facet] = makeFacet(facet, mask);
108 0 : AlwaysAssert(!facetMaskImages_p[facet].null(), AipsError);
109 0 : AlwaysAssert(MFCleanImageSkyModel::addMask(facet, *facetMaskImages_p[facet]),
110 : AipsError);
111 : }
112 0 : return true;
113 : }
114 : else {
115 : // All other images are added as is
116 0 : return MFCleanImageSkyModel::addMask(image+nfacets_p-1, mask);
117 : }
118 : };
119 :
120 0 : Bool WFCleanImageSkyModel::addResidual(Int image,
121 : ImageInterface<Float>& residual)
122 : {
123 : // Create the facet images and copy the relevant region from the
124 : // original image
125 0 : AlwaysAssert(image>-1, AipsError);
126 0 : if(residualImage_p.null()) {
127 0 : residualImage_p=CountedPtr<ImageInterface<Float> >(&residual, false);
128 0 : facetResidualImages_p.resize(nfacets_p);
129 0 : for (Int facet=0;facet<nfacets_p;facet++) {
130 0 : facetResidualImages_p[facet] = makeFacet(facet, residual);
131 0 : AlwaysAssert(!facetResidualImages_p[facet].null(), AipsError);
132 0 : AlwaysAssert(MFCleanImageSkyModel::addResidual(facet, *facetResidualImages_p[facet]),
133 : AipsError);
134 : }
135 0 : return true;
136 : }
137 : else {
138 : // All other images are added as is
139 0 : return MFCleanImageSkyModel::addResidual(image+nfacets_p-1, residual);
140 : }
141 : };
142 :
143 0 : WFCleanImageSkyModel::~WFCleanImageSkyModel()
144 : {
145 : // Now destroy all the facet images
146 : /*
147 : if(imageImage_p !=0){
148 : for (Int facet=0;facet<nfacets_p;facet++) {
149 : if(facetImages_p[facet]); delete facetImages_p[facet];
150 : facetImages_p[facet]=0;
151 : if(Int(facetMaskImages_p.nelements())>facet) {
152 : if(facetMaskImages_p[facet]) delete facetMaskImages_p[facet];
153 : facetMaskImages_p[facet]=0;
154 : }
155 : if(Int(facetResidualImages_p.nelements())>facet) {
156 : if(facetResidualImages_p[facet]) delete facetResidualImages_p[facet];
157 : facetResidualImages_p[facet]=0;
158 : }
159 : }
160 : }
161 : */
162 0 : };
163 :
164 : // Clean solver: just call the MF solver
165 0 : Bool WFCleanImageSkyModel::solve(SkyEquation& se) {
166 :
167 0 : LogIO os(LogOrigin("WFCleanImageSkyModel","solve"));
168 :
169 : os << "Starting wide-field clean with " << nfacets_p << " facets"
170 0 : << " and " << this->numberOfModels()-nfacets_p << " outliers" << LogIO::POST;
171 0 : maskImage_p=0; // the masks will be have to be renewed when continuing
172 : // The MFCleanImageSkyModel solve cleans all images in
173 : // parallel.
174 0 : if(!MFCleanImageSkyModel::solve(se)) {
175 0 : os << "Wide-field clean apparently failed to reach threshold" << LogIO::POST;
176 0 : return false;
177 : }
178 :
179 0 : return(true);
180 0 : };
181 :
182 0 : ImageInterface<Float>& WFCleanImageSkyModel::getResidual(Int trueMod){
183 :
184 0 : if(trueMod==0){
185 0 : if(residualImage_p.null()){
186 0 : residualImage_p=new TempImage<Float>(imageImage_p->shape(), imageImage_p->coordinates());
187 0 : facetResidualImages_p.resize(nfacets_p);
188 0 : for (Int facet=0;facet<nfacets_p; ++facet) {
189 0 : facetResidualImages_p[facet] = makeFacet(facet, *residualImage_p);
190 0 : AlwaysAssert(!facetResidualImages_p[facet].null(), AipsError);
191 : }
192 :
193 : }
194 0 : for (Int facet=0;facet<nfacets_p;++facet)
195 0 : facetResidualImages_p[facet]->copyData(residual(facet));
196 0 : return *residualImage_p;
197 :
198 : }
199 : else{
200 : //return the outlying Multi field residuals
201 0 : return residual(trueMod+nfacets_p-1);
202 :
203 : }
204 :
205 :
206 :
207 : }
208 :
209 : SubImage<Float>*
210 0 : WFCleanImageSkyModel::makeFacet(Int facet, ImageInterface<Float>& image)
211 : {
212 :
213 0 : LogIO os(LogOrigin("WFCleanImageSkyModel","makeFacet"));
214 :
215 :
216 : // Make the output image
217 0 : Slicer imageSlicer;
218 :
219 : // Have to fiddle about with the shape (double it in x and y)
220 : // and the coordinates (change the ref pixel).
221 0 : IPosition facetShape(image.shape());
222 0 : CoordinateSystem coords=image.coordinates();
223 0 : Int directionIndex=coords.findCoordinate(Coordinate::DIRECTION);
224 0 : AlwaysAssert(directionIndex>=0, AipsError);
225 : DirectionCoordinate
226 0 : facetDirCoord=coords.directionCoordinate(directionIndex);
227 0 : if(!makeSlicers(facet, image.shape(), facetShape, imageSlicer)) {
228 0 : os << LogIO::SEVERE << "Cannot create facet image" << LogIO::EXCEPTION;
229 : }
230 : // Now we have all the information and we can create the facet image
231 0 : SubImage<Float>* facetImage = new SubImage<Float>(image, imageSlicer, true);
232 0 : facetImage->setMiscInfo(image.miscInfo());
233 0 : facetImage->setUnits(image.units());
234 :
235 0 : return facetImage;
236 0 : }
237 :
238 : Bool
239 0 : WFCleanImageSkyModel::makeSlicers(const Int facet, const IPosition& imageShape,
240 : IPosition& facetShape,
241 : Slicer& imageSlicer)
242 : {
243 0 : LogIO os(LogOrigin("WFCleanImageSkyModel","makeSlicers"));
244 :
245 0 : if((facet>(nfacets_p-1))||(facet<0)) {
246 0 : os << LogIO::SEVERE << "Illegal facet " << facet << LogIO::POST;
247 0 : return false;
248 : }
249 :
250 0 : IPosition imageBlc(imageShape.nelements(), 0);
251 0 : IPosition imageTrc(imageShape.nelements(), 0);
252 0 : IPosition inc=IPosition(imageShape.nelements(), 1);
253 :
254 0 : facetShape=imageShape;
255 0 : IPosition facetBlc(facetShape.nelements(), 0);
256 0 : IPosition facetTrc(facetShape.nelements(), 0);
257 :
258 0 : Int facetx = facet % facets_p;
259 0 : Int facety = (facet - facetx) / facets_p;
260 0 : Int sizex = imageShape(0) / facets_p;
261 0 : Int sizey = imageShape(1) / facets_p;
262 :
263 : // Make the image Slicer
264 0 : imageBlc(0) = facetx * sizex; imageTrc(0) = imageBlc(0) + sizex - 1;
265 0 : imageBlc(1) = facety * sizey; imageTrc(1) = imageBlc(1) + sizey - 1;
266 : {
267 0 : for (uInt i=2;i<imageShape.nelements();i++) {
268 0 : imageTrc(i)=imageBlc(i)+imageShape(i)-1;
269 0 : if(imageTrc(i) > (imageShape(i)-1)) imageTrc(i) = imageShape(i) -1;
270 : }
271 : }
272 :
273 0 : LCBox::verify(imageBlc, imageTrc, inc, imageShape);
274 :
275 0 : imageSlicer=Slicer (imageBlc, imageTrc, inc, Slicer::endIsLast);
276 :
277 : // Make the facet Slicer to place the imageBlc at the
278 : // end of first quarter of the facet and the imageTrc at the
279 : // start of the last quarter.
280 0 : facetShape(0) = sizex;
281 0 : facetBlc(0) = 0; facetTrc(0) = facetBlc(0) + sizex - 1;
282 0 : facetShape(1) = sizey;
283 0 : facetBlc(1) = 0; facetTrc(1) = facetBlc(1) + sizey - 1;
284 : {
285 0 : for (uInt i=2;i<facetShape.nelements();i++) {
286 0 : facetTrc(i)=facetBlc(i)+facetShape(i)-1;
287 0 : if(facetTrc(i) > (facetShape(i)-1)) facetTrc(i) = facetShape(i) -1;
288 : }
289 : }
290 :
291 0 : LCBox::verify(facetBlc, facetTrc, inc, facetShape);
292 :
293 : os << LogIO::DEBUGGING << "Facet " << facet+1
294 0 : << " : from " << imageBlc+1<< " to " << imageTrc+1 << LogIO::POST;
295 :
296 0 : return true;
297 0 : }
298 :
299 :
300 :
301 : } //# NAMESPACE CASA - END
302 :
303 :
304 :
|
