Line data Source code
1 : /*
2 : * LinearMosaic.cc
3 : //# Copyright (C) 2015
4 : //# Associated Universities, Inc. Washington DC, USA.
5 : //#
6 : //# This program is free software; you can redistribute it and/or modify it
7 : //# under the terms of the GNU General Public License as published by the Free
8 : //# Software Foundation; either version 2 of the License, or (at your option)
9 : //# any later version.
10 : //#
11 : //# This program is distributed in the hope that it will be useful, but WITHOUT
12 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
13 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
14 : //# more details.
15 : //#
16 : //# You should have received a copy of the GNU General Public License along
17 : //# with this program; if not, write to the Free Software Foundation, Inc.,
18 : //# 675 Massachusetts Ave, Cambridge, MA 02139, USA.
19 : //#
20 : //# Correspondence concerning AIPS++ should be addressed as follows:
21 : //# Internet email: casa-feedback@nrao.edu.
22 : //# Postal address: AIPS++ Project Office
23 : //# National Radio Astronomy Observatory
24 : //# 520 Edgemont Road
25 : //# Charlottesville, VA 22903-2475 USA
26 : * Created on: Feb 24, 2015
27 : * Author: kgolap
28 : */
29 :
30 : #include <casacore/coordinates/Coordinates/CoordinateUtil.h>
31 : #include <casacore/coordinates/Coordinates/CoordinateSystem.h>
32 : #include <casacore/coordinates/Coordinates/DirectionCoordinate.h>
33 : #include <casacore/images/Images/TempImage.h>
34 : #include <casacore/images/Images/PagedImage.h>
35 : #include <casacore/images/Images/ImageRegrid.h>
36 : #include <casacore/images/Images/SubImage.h>
37 : #include <casacore/images/Regions/WCBox.h>
38 : #include <synthesis/MeasurementEquations/Imager.h>
39 : #include <synthesis/MeasurementEquations/LinearMosaic.h>
40 : using namespace std;
41 :
42 : using namespace casacore;
43 : namespace casa { //# NAMESPACE CASA - BEGIN
44 :
45 :
46 0 : LinearMosaic::LinearMosaic (): outImage_p(NULL), outWgt_p(NULL), outImName_p(""), outWgtName_p(""), weightType_p(1), linmosType_p(2){
47 :
48 :
49 0 : }
50 0 : LinearMosaic::LinearMosaic(const String outim, const String outwgt, const MDirection& imcen, const Int nx,
51 0 : const Int ny, const Quantity cellx, const Quantity celly, const Int linmostype) :outImage_p(NULL), outWgt_p(NULL), linmosType_p(linmostype) {
52 0 : outImName_p=outim;
53 : ///Null image so assigning same weight type as requested
54 0 : imageWeightType_p=linmostype;
55 0 : weightType_p=linmostype;
56 :
57 0 : outWgtName_p= outwgt=="" ? outim+String(".weight") : outwgt;
58 0 : nx_p=nx;
59 0 : ny_p=ny;
60 0 : imcen_p=imcen;
61 0 : cellx_p=cellx;
62 0 : celly_p=celly;
63 :
64 0 : }
65 0 : LinearMosaic::LinearMosaic(const String outim, const String outwgt, const MDirection& imcen, const Int nx, const Int ny,
66 : Vector<CountedPtr<ImageInterface<Float> > >& ims,
67 0 : Vector<CountedPtr<ImageInterface<Float> > >& wgtims, const Int linmostype) : linmosType_p(linmostype){
68 0 : Int nchan=ims[0]->shape()[3];
69 0 : Int npol=ims[0]->shape()[2];
70 0 : for (uInt k=0; k< ims.nelements(); ++k){
71 0 : if(nchan != ims[k]->shape()[3] || nchan != wgtims[k]->shape()[3])
72 0 : ThrowCc("images should have the same number of channels");
73 0 : if(npol != ims[k]->shape()[2] || npol != (wgtims[k])->shape()[2])
74 0 : ThrowCc("images should have the same number of polarization planes");
75 : }
76 0 : CoordinateSystem cs=ims[0]->coordinates();
77 0 : makeEmptyImage(outim, cs, imcen, nx, ny,npol, nchan);
78 0 : makeEmptyImage(outwgt, cs, imcen, nx, ny,npol, nchan);
79 0 : PagedImage<Float> outdiskim(outim);
80 0 : PagedImage<Float> outdiskwgt(outwgt);
81 0 : setOutImages(outdiskim, outdiskwgt, 2, 2);
82 0 : makeMosaic(ims, wgtims);
83 :
84 0 : }
85 0 : void LinearMosaic::setlinmostype(const Int linmostype){
86 0 : linmosType_p=linmostype;
87 0 : }
88 0 : Bool LinearMosaic::makeMosaic(ImageInterface<Float>& outim, ImageInterface<Float>& outwgt,
89 : Vector<CountedPtr<ImageInterface<Float> > >& ims,
90 : Vector<CountedPtr<ImageInterface<Float> > >& wgtims){
91 :
92 0 : Bool retval=true;
93 0 : if(ims.nelements() != wgtims.nelements())
94 0 : ThrowCc("Unequal number of images and weight images ");
95 0 : for (uInt k=0; k < ims.nelements(); ++k){
96 0 : retval=retval && addOnToImage(outim, outwgt, *(ims[k]), *(wgtims[k]), (k==ims.nelements()-1));
97 : }
98 0 : return retval;
99 :
100 : }
101 0 : Bool LinearMosaic::makeMosaic(
102 : Vector<CountedPtr<ImageInterface<Float> > >& ims,
103 : Vector<CountedPtr<ImageInterface<Float> > >& wgtims){
104 :
105 0 : Bool retval=true;
106 0 : if(outImage_p.null() && outImName_p=="")
107 0 : ThrowCc("No output image or weight image defined");
108 0 : if(outImage_p.null())
109 0 : createOutImages(ims[0]->coordinates(), ims[0]->shape()[2], ims[0]->shape()[3] );
110 0 : if(ims.nelements() != wgtims.nelements())
111 0 : ThrowCc("Unequal number of images and weight images ");
112 0 : for (uInt k=0; k < ims.nelements(); ++k){
113 0 : retval=retval && addOnToImage(*outImage_p, *outWgt_p, *(ims[k]), *(wgtims[k]), (k==ims.nelements()-1));
114 0 : ims[k]=0;
115 0 : wgtims[k]=0;
116 : }
117 0 : return retval;
118 :
119 : }
120 0 : void LinearMosaic::setOutImages(ImageInterface<Float>& outim, ImageInterface<Float>& outwgt, const Int imageWeightType, const Int weightType){
121 0 : outImage_p=CountedPtr<ImageInterface<Float> >(&outim, false);
122 0 : outWgt_p = CountedPtr<ImageInterface<Float> >(&outwgt, false);
123 0 : imageWeightType_p=imageWeightType;
124 0 : weightType_p=weightType;
125 0 : }
126 0 : void LinearMosaic::setOutImages(const String& outim, const String& outwgt, const Int imageWeightType, const Int weightType){
127 0 : outImage_p=new PagedImage<Float>(outim);
128 0 : outWgt_p = new PagedImage<Float>(outwgt);
129 0 : imageWeightType_p=imageWeightType;
130 0 : weightType_p=weightType;
131 0 : }
132 0 : Bool LinearMosaic::addOnToImage(ImageInterface<Float>& outim, ImageInterface<Float>& outwgt, const ImageInterface<Float>& inIm,
133 : const ImageInterface<Float>& inWgt, Bool unWeightOutIm){
134 :
135 0 : if(inIm.shape() != inWgt.shape())
136 0 : ThrowCc("Image "+ inIm.name() + " and Weight image "+inWgt.name() +". have to be similar...please regrid appropriately");
137 0 : Double meminMB=Double(HostInfo::memoryTotal(true))/1024.0;
138 :
139 : //if(!outImIsWeighted){
140 : {
141 : // cerr << "limosType " << linmosType_p << " imageWeightType " << imageWeightType_p << " weightType " << weightType_p << endl;
142 0 : if(linmosType_p==2){
143 0 : if(weightType_p==1)
144 0 : outwgt.copyData((LatticeExpr<Float>)(outwgt*outwgt));
145 0 : weightType_p=2;
146 0 : if(imageWeightType_p==1)
147 0 : outim.copyData((LatticeExpr<Float>)(outim*sqrt(outwgt)));
148 0 : else if( imageWeightType_p==0)
149 0 : outim.copyData((LatticeExpr<Float>)(outim*outwgt));
150 0 : imageWeightType_p=2;
151 :
152 : }
153 : //cerr << " imageWeightType " << imageWeightType_p << " weightType " << weightType_p << endl;
154 0 : if(linmosType_p==1){
155 0 : if(weightType_p==2)
156 0 : outwgt.copyData((LatticeExpr<Float>)(sqrt(abs(outwgt))));
157 0 : weightType_p=1;
158 0 : if(imageWeightType_p==0)
159 0 : outim.copyData((LatticeExpr<Float>)(outim*outwgt));
160 0 : else if( imageWeightType_p==2)
161 0 : outim.copyData((LatticeExpr<Float>)(iif(outwgt > (0.0),
162 0 : (outim/outwgt), 0)));
163 0 : imageWeightType_p=1;
164 :
165 : }
166 : }
167 :
168 0 : CoordinateSystem incs=inIm.coordinates();
169 :
170 :
171 :
172 0 : IPosition iblc(inIm.shape().nelements(),0);
173 0 : IPosition itrc=inIm.shape();
174 : // cerr << "itrc " << itrc << endl;
175 0 : itrc -=1;
176 :
177 :
178 :
179 : //cerr << "blc " << iblc << " trc " << itrc << endl;
180 0 : LCBox lbox(iblc, itrc, inIm.shape());
181 0 : WCBox wbox(lbox, incs);
182 : //cerr << "wbox " << wbox.toRecord("") << endl;
183 0 : ImageRegion imagreg(wbox );
184 0 : SubImage<Float> subOutIm;
185 0 : SubImage<Float> subOutWgt;
186 : try{
187 :
188 0 : subOutIm=SubImage<Float>(outim, imagreg, true);
189 0 : subOutWgt=SubImage<Float>(outwgt, imagreg, true);
190 : }
191 0 : catch(...){
192 : //Failed to make a subimage let us use the full image
193 : //cerr << "Failed to make subImage " << x.what()<< endl;
194 0 : subOutIm=SubImage<Float>(outim, true);
195 0 : subOutWgt=SubImage<Float>(outwgt, true);
196 :
197 0 : }
198 0 : TempImage<Float> fullImage(subOutWgt.shape(), subOutIm.coordinates(), meminMB/8.0);
199 0 : TempImage<Float> fullWeight(subOutWgt.shape(), subOutIm.coordinates(), meminMB/8.0);
200 0 : fullImage.set(0.0);
201 0 : fullWeight.set(0.0);
202 0 : ImageRegrid<Float> regridder;
203 : {
204 :
205 0 : if(linmosType_p==2){
206 0 : TempImage<Float> trueWeightIm(inWgt.shape(), inWgt.coordinates(), meminMB/8.0);
207 0 : if(inWgt.getDefaultMask() != ""){
208 0 : Imager::copyMask(trueWeightIm, inWgt, inWgt.getDefaultMask());
209 0 : fullWeight.makeMask(inWgt.getDefaultMask(), true, true, true, true);
210 : }
211 0 : trueWeightIm.copyData((LatticeExpr<Float>)(inWgt*inWgt));
212 :
213 0 : regridder.regrid( fullWeight, Interpolate2D::LINEAR,
214 0 : IPosition(2,0,1), trueWeightIm);
215 0 : TempImage<Float> inWeightedIm(inIm.shape(), inIm.coordinates(), meminMB/8.0);
216 0 : if(inIm.getDefaultMask() != ""){
217 0 : Imager::copyMask(inWeightedIm, inIm, inIm.getDefaultMask());
218 0 : fullImage.makeMask(inIm.getDefaultMask(), true, true, true, true);
219 : }
220 :
221 0 : inWeightedIm.copyData((LatticeExpr<Float>)(inIm*inWgt));
222 0 : ImageRegrid<Float> regridder;
223 0 : regridder.regrid( fullImage, Interpolate2D::LINEAR,
224 0 : IPosition(2,0,1), inWeightedIm);
225 :
226 0 : }
227 0 : else if (linmosType_p==1){
228 0 : if(inIm.getDefaultMask() != ""){
229 0 : fullImage.makeMask(inIm.getDefaultMask(), true, true, true, true);
230 : }
231 0 : if(inWgt.getDefaultMask() != ""){
232 0 : fullWeight.makeMask(inWgt.getDefaultMask(), true, true, true, true);
233 : }
234 0 : regridder.regrid( fullWeight, Interpolate2D::LINEAR,
235 0 : IPosition(2,0,1), inWgt);
236 0 : regridder.regrid( fullImage, Interpolate2D::LINEAR,
237 0 : IPosition(2,0,1), inIm);
238 : }
239 : }
240 :
241 0 : subOutWgt.copyData((LatticeExpr<Float>)(subOutWgt+fullWeight));
242 :
243 0 : subOutIm.copyData((LatticeExpr<Float>)(subOutIm+fullImage));
244 : //if(wMax > 0.0){
245 : // outim.copyData((LatticeExpr<Float>)(outim/wMax));
246 : // outwgt.copyData((LatticeExpr<Float>)(outwgt/wMax));
247 : //}
248 0 : if(unWeightOutIm){
249 : //LatticeExprNode outmax0 = max( outim );
250 :
251 : //cerr << "Bef max of out " << outmax0.getFloat() << endl;
252 :
253 0 : LatticeExprNode elmax = max( outwgt );
254 :
255 : //Cannot trust anything below 1e-6
256 0 : Float wMinLimit = 1e-6* elmax.getFloat();
257 0 : outim.copyData((LatticeExpr<Float>)(iif(outwgt > (wMinLimit),
258 0 : (outim/outwgt), 0)));
259 : //LatticeExprNode outmax = max( outim );
260 :
261 : //cerr << "max of out " << outmax.getFloat() << endl;
262 0 : imageWeightType_p=0;
263 0 : }
264 :
265 0 : return true;
266 0 : }
267 :
268 :
269 0 : void LinearMosaic::saultWeightImage(const String& outimname, const Float& fracPeakWgt){
270 0 : if(outImage_p.null())
271 0 : ThrowCc("Mosaic image and weight must be set");
272 0 : PagedImage<Float> outdiskim(outImage_p->shape(), outImage_p->coordinates(), outimname);
273 0 : LatticeExprNode elmax = max( *outWgt_p );
274 0 : Float wMax = elmax.getFloat();
275 0 : wMax *=fracPeakWgt;
276 0 : LatticeExpr<Float> weightMath;
277 0 : if(imageWeightType_p==0){
278 0 : if(weightType_p==2)
279 0 : weightMath=(LatticeExpr<Float>)(iif((*outWgt_p) > (wMax),
280 0 : (*outImage_p), (*outImage_p)*(*outWgt_p)));
281 0 : else if(weightType_p==1)
282 0 : weightMath=(LatticeExpr<Float>)(iif((*outWgt_p) > (wMax),
283 0 : (*outImage_p), (*outImage_p)*(*outWgt_p)*(*outWgt_p)));
284 : }
285 0 : if(imageWeightType_p==1){
286 0 : if( weightType_p==2)
287 0 : weightMath=(LatticeExpr<Float>)(iif((*outWgt_p) > (wMax),
288 0 : (*outImage_p)/sqrt(*outWgt_p), (*outImage_p)*sqrt(*outWgt_p)));
289 0 : else if(weightType_p==1)
290 0 : weightMath=(LatticeExpr<Float>)(iif((*outWgt_p) > (wMax),
291 0 : (*outImage_p)/(*outWgt_p), (*outImage_p)*(*outWgt_p)));
292 : }
293 0 : if(imageWeightType_p==2){
294 0 : if( weightType_p==2)
295 0 : weightMath=(LatticeExpr<Float>)(iif((*outWgt_p) > (wMax),
296 0 : (*outImage_p)/(*outWgt_p), (*outImage_p)));
297 0 : else if(weightType_p==1)
298 0 : weightMath=(LatticeExpr<Float>)(iif((*outWgt_p) > (wMax),
299 0 : (*outImage_p)/((*outWgt_p)*(*outWgt_p)), (*outImage_p)));
300 : }
301 :
302 0 : outdiskim.copyData(weightMath);
303 :
304 :
305 0 : }
306 0 : void LinearMosaic::createOutImages(const CoordinateSystem& cs, const Int npol, const Int nchan ){
307 0 : makeEmptyImage(outImName_p, cs, imcen_p, nx_p, ny_p, npol, nchan );
308 0 : makeEmptyImage(outWgtName_p, cs, imcen_p, nx_p, ny_p, npol, nchan );
309 0 : outImage_p=new PagedImage<Float>(outImName_p);
310 0 : outWgt_p= new PagedImage<Float>(outWgtName_p);
311 0 : }
312 0 : void LinearMosaic::makeEmptyImage(const String imagename, const CoordinateSystem& cs, const MDirection& imcen, const Int nx, const Int ny,
313 : const Int npol, const Int nchan){
314 0 : CoordinateSystem outcs=cs;
315 0 : Int dirAx=outcs.findCoordinate(Coordinate::DIRECTION);
316 0 : DirectionCoordinate dc=outcs.directionCoordinate(dirAx);
317 0 : Vector<Double> incr=dc.increment();
318 0 : String elunit=dc.worldAxisUnits()[0];
319 0 : if(cellx_p.getValue() != 0.0)
320 0 : incr[0]=-cellx_p.getValue(elunit);
321 0 : if(celly_p.getValue() != 0.0)
322 0 : incr[1]=celly_p.getValue(elunit);
323 0 : dc.setIncrement(incr);
324 0 : Vector<Double> cenVec(2);
325 0 : cenVec[0]=Double(nx)/2.0; cenVec[1]=Double(ny)/2.0;
326 0 : dc.setReferencePixel(cenVec);
327 : MDirection::Types eltype;
328 0 : MDirection::getType(eltype, imcen.getRefString());
329 0 : dc.setReferenceValue(imcen.getAngle().getValue(elunit));
330 0 : dc.setReferenceFrame(eltype);
331 0 : outcs.replaceCoordinate(dc, dirAx);
332 0 : PagedImage<Float> outdiskim(IPosition(4, nx, ny, npol, nchan), outcs, imagename);
333 0 : outdiskim.set(0.0);
334 0 : }
335 :
336 : using namespace casacore;
337 : } //end namespace casa
|
