Line data Source code
1 : //# Feather.cc: Implementation of Feather.h
2 : //# Copyright (C) 1997-2013
3 : //# Associated Universities, Inc. Washington DC, USA.
4 : //#
5 : //# This program is free software; you can redistribute it and/or modify it
6 : //# under the terms of the GNU General Public License as published by the Free
7 : //# Software Foundation; either version 2 of the License, or (at your option)
8 : //# any later version.
9 : //#
10 : //# This program 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 General Public License for
13 : //# more details.
14 : //#
15 : //# You should have received a copy of the GNU General Public License along
16 : //# with this program; if not, write to the Free Software Foundation, Inc.,
17 : //# 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 : //# $kgolap$
27 : #include <cmath>
28 : #include <synthesis/MeasurementEquations/Feather.h>
29 : #include <synthesis/MeasurementEquations/Imager.h>
30 : #include <synthesis/TransformMachines/StokesImageUtil.h>
31 : #include <casacore/casa/OS/HostInfo.h>
32 :
33 : #include <casacore/casa/Arrays/Matrix.h>
34 : #include <casacore/casa/Arrays/ArrayMath.h>
35 : #include <casacore/casa/Arrays/ArrayLogical.h>
36 : #include <casacore/casa/Arrays/IPosition.h>
37 : #include <iostream>
38 : #include <casacore/casa/Logging.h>
39 : #include <casacore/casa/Logging/LogIO.h>
40 : #include <casacore/casa/Logging/LogMessage.h>
41 : #include <casacore/casa/Logging/LogSink.h>
42 : #include <casacore/scimath/Mathematics/MathFunc.h>
43 : #include <casacore/tables/TaQL/ExprNode.h>
44 : #include <casacore/casa/Utilities/Assert.h>
45 : #include <casacore/casa/Arrays/ArrayMath.h>
46 : #include <casacore/casa/Arrays/Slice.h>
47 : #include <casacore/images/Images/TempImage.h>
48 : #include <casacore/images/Images/ImageInterface.h>
49 : #include <casacore/images/Images/PagedImage.h>
50 : #include <casacore/images/Images/ImageRegrid.h>
51 : #include <casacore/images/Images/ImageUtilities.h>
52 : #include <synthesis/TransformMachines/PBMath.h>
53 : #include <casacore/lattices/LEL/LatticeExpr.h>
54 : #include <casacore/lattices/LatticeMath/LatticeFFT.h>
55 : #include <casacore/coordinates/Coordinates/CoordinateSystem.h>
56 : #include <casacore/coordinates/Coordinates/DirectionCoordinate.h>
57 : #include <casacore/coordinates/Coordinates/SpectralCoordinate.h>
58 : #include <casacore/coordinates/Coordinates/StokesCoordinate.h>
59 : #include <casacore/coordinates/Coordinates/CoordinateUtil.h>
60 : #include <casacore/coordinates/Coordinates/Projection.h>
61 : #include <casacore/coordinates/Coordinates/ObsInfo.h>
62 :
63 : #include <components/ComponentModels/GaussianDeconvolver.h>
64 :
65 : using namespace casacore;
66 : namespace casa { //# NAMESPACE CASA - BEGIN
67 :
68 :
69 113 : Feather::Feather(): dishDiam_p(-1.0), cweightCalced_p(false), cweightApplied_p(false), sdScale_p(1.0){
70 113 : highIm_p=NULL;
71 113 : lowIm_p=NULL;
72 113 : lowImOrig_p=NULL;
73 113 : cwImage_p=NULL;
74 113 : cwHighIm_p=NULL;
75 113 : }
76 0 : Feather::Feather(const ImageInterface<Float>& SDImage, const ImageInterface<Float>& INTImage, Float sdscale) : dishDiam_p(-1.0), cweightCalced_p(false), cweightApplied_p(false), sdScale_p(sdscale){
77 :
78 0 : setINTImage(INTImage);
79 0 : lowImOrig_p=NULL;
80 0 : setSDImage(SDImage);
81 :
82 0 : cwImage_p=NULL;
83 0 : cwHighIm_p=NULL;
84 :
85 0 : }
86 113 : Feather::~Feather(){
87 113 : highIm_p=NULL;
88 113 : lowIm_p=NULL;
89 113 : cwImage_p=NULL;
90 113 : }
91 :
92 113 : void Feather::setSDImage(const ImageInterface<Float>& SDImage){
93 339 : LogIO os(LogOrigin("Feather", "setSDImage()", WHERE));
94 113 : if(highIm_p.null())
95 0 : throw(AipsError("High res image has to be defined before SD image"));
96 226 : ImageInfo lowInfo=SDImage.imageInfo();
97 113 : lBeam_p=lowInfo.restoringBeam();
98 113 : if(lBeam_p.isNull())
99 0 : throw(AipsError("No Single dish restoring beam info in image"));
100 226 : CoordinateSystem csyslow=SDImage.coordinates();
101 226 : CountedPtr<ImageInterface<Float> > sdcopy;
102 113 : sdcopy=new TempImage<Float>(SDImage.shape(), csyslow);
103 113 : (*sdcopy).copyData(SDImage);
104 113 : ImageUtilities::copyMiscellaneous(*sdcopy, SDImage);
105 113 : if(SDImage.getDefaultMask() != "")
106 33 : Imager::copyMask(*sdcopy, SDImage, SDImage.getDefaultMask());
107 113 : std::unique_ptr<ImageInterface<Float> > copyPtr;
108 113 : std::unique_ptr<ImageInterface<Float> > copyPtr2;
109 :
110 :
111 113 : Vector<Stokes::StokesTypes> stokesvec;
112 113 : if(CoordinateUtil::findStokesAxis(stokesvec, csyslow) <0){
113 0 : CoordinateUtil::addIAxis(csyslow);
114 :
115 0 : ImageUtilities::addDegenerateAxes (os, copyPtr, *sdcopy, "",
116 : false, false,"I", false, false,
117 : true);
118 0 : sdcopy=CountedPtr<ImageInterface<Float> >(copyPtr.get(), false);
119 :
120 : }
121 113 : if(CoordinateUtil::findSpectralAxis(csyslow) <0){
122 0 : CoordinateUtil::addFreqAxis(csyslow);
123 0 : ImageUtilities::addDegenerateAxes (os, copyPtr2, *sdcopy, "",
124 : false, true,
125 : "", false, false,
126 : true);
127 0 : sdcopy=CountedPtr<ImageInterface<Float> >(copyPtr2.get(), false);
128 : }
129 113 : lowIm_p=new TempImage<Float>(highIm_p->shape(), csysHigh_p);
130 : // regrid the single dish image
131 : {
132 226 : Vector<Int> dirAxes=CoordinateUtil::findDirectionAxes(csysHigh_p);
133 226 : IPosition axes(3,dirAxes(0),dirAxes(1),2);
134 113 : Int spectralAxisIndex=CoordinateUtil::findSpectralAxis(csysHigh_p);
135 113 : axes(2)=spectralAxisIndex;
136 113 : if(sdcopy->getDefaultMask() != "")
137 33 : lowIm_p->makeMask(sdcopy->getDefaultMask(), true, true, true, true);
138 :
139 226 : ImageRegrid<Float> ir;
140 113 : ir.regrid(*lowIm_p, Interpolate2D::LINEAR, axes, *sdcopy);
141 : }
142 : /*if(sdcopy->getDefaultMask() != ""){
143 : //Imager::copyMask(*lowIm_p, *sdcopy, sdcopy->getDefaultMask());
144 : lowIm_p->makeMask(sdcopy->getDefaultMask(), true, true);
145 : ImageUtilities::copyMask(*lowIm_p, *sdcopy,sdcopy->getDefaultMask() , sdcopy->getDefaultMask(), AxesSpecifier());
146 : lowIm_p->setDefaultMask(sdcopy->getDefaultMask());
147 :
148 : }
149 : */
150 :
151 113 : if(lowImOrig_p.null()){
152 113 : lowImOrig_p=new TempImage<Float>(lowIm_p->shape(), lowIm_p->coordinates(),0);
153 113 : lowImOrig_p->copyData(*lowIm_p);
154 113 : lBeamOrig_p=lBeam_p;
155 : }
156 113 : cweightCalced_p=false;
157 113 : }
158 :
159 0 : void Feather::convolveINT(const GaussianBeam& newHighBeam){
160 0 : GaussianBeam toBeUsed(Quantity(0.0, "arcsec"),Quantity(0.0, "arcsec"), Quantity(0.0, "deg")) ;
161 0 : Bool retval=true;
162 : try {
163 : //cerr << "highBeam " << hBeam_p.getMajor() << " " << hBeam_p.getMinor() << " " << hBeam_p.getPA() << endl;
164 : retval=
165 0 : GaussianDeconvolver::deconvolve(toBeUsed, newHighBeam, hBeam_p);
166 : //cerr << "beam to be used " << toBeUsed.getMajor() << " " << toBeUsed.getMinor() << " " << toBeUsed.getPA() << endl;
167 : }
168 0 : catch (const AipsError& x) {
169 0 : if(toBeUsed.getMajor("arcsec")==0.0)
170 0 : throw(AipsError("new Beam may be smaller than the beam of original Interferometer image"));
171 : }
172 : try{
173 0 : StokesImageUtil::Convolve(*highIm_p, toBeUsed, true);
174 : }
175 0 : catch(const AipsError& x){
176 0 : throw(AipsError("Could not convolve INT image for some reason; try a lower resolution may be"));
177 : }
178 0 : hBeam_p=newHighBeam;
179 :
180 : //need to redo feather application
181 0 : cweightApplied_p=false;
182 : (void)retval; // avoid compiler warning
183 0 : }
184 :
185 113 : void Feather::setINTImage(const ImageInterface<Float>& INTImage){
186 339 : LogIO os(LogOrigin("Feather", "setINTImage()", WHERE));
187 226 : ImageInfo highInfo=INTImage.imageInfo();
188 113 : hBeam_p=highInfo.restoringBeam();
189 113 : if(hBeam_p.isNull())
190 0 : throw(AipsError("No high-resolution restoring beam info in image"));
191 113 : csysHigh_p=INTImage.coordinates();
192 339 : CountedPtr<ImageInterface<Float> > intcopy=new TempImage<Float>(INTImage.shape(), csysHigh_p);
193 113 : (*intcopy).copyData(INTImage);
194 113 : ImageUtilities::copyMiscellaneous(*intcopy, INTImage);
195 113 : if(INTImage.getDefaultMask() != "")
196 16 : Imager::copyMask(*intcopy, INTImage, INTImage.getDefaultMask());
197 113 : std::unique_ptr<ImageInterface<Float> > copyPtr;
198 113 : std::unique_ptr<ImageInterface<Float> > copyPtr2;
199 :
200 :
201 226 : Vector<Stokes::StokesTypes> stokesvec;
202 113 : if(CoordinateUtil::findStokesAxis(stokesvec, csysHigh_p) <0){
203 0 : CoordinateUtil::addIAxis(csysHigh_p);
204 0 : ImageUtilities::addDegenerateAxes (os, copyPtr, *intcopy, "",
205 : false, false,"I", false, false,
206 : true);
207 0 : intcopy=CountedPtr<ImageInterface<Float> >(copyPtr.get(), false);
208 :
209 : }
210 113 : if(CoordinateUtil::findSpectralAxis(csysHigh_p) <0){
211 0 : CoordinateUtil::addFreqAxis(csysHigh_p);
212 0 : ImageUtilities::addDegenerateAxes (os, copyPtr2, *intcopy, "",
213 : false, true,
214 : "", false, false,
215 : true);
216 0 : intcopy=CountedPtr<ImageInterface<Float> >(copyPtr2.get(), false);
217 : }
218 :
219 :
220 113 : highIm_p=new TempImage<Float>(intcopy->shape(), csysHigh_p);
221 :
222 113 : highIm_p->copyData(*intcopy);
223 113 : ImageUtilities::copyMiscellaneous(*highIm_p, *intcopy);
224 113 : String maskname=intcopy->getDefaultMask();
225 113 : if(maskname != ""){
226 16 : Imager::copyMask(*highIm_p, *intcopy, maskname);
227 :
228 : }
229 113 : cweightCalced_p=false;
230 :
231 113 : }
232 :
233 106 : Bool Feather::setEffectiveDishDiam(const Float xdiam, const Float ydiam){
234 106 : dishDiam_p=ydiam >0 ? min(xdiam, ydiam) : xdiam;
235 : {//reset to the original SD image
236 106 : lowIm_p->copyData(*lowImOrig_p);
237 106 : lBeam_p=lBeamOrig_p;
238 : }
239 : //Change the beam of SD image
240 106 : Double freq = worldFreq(lowIm_p->coordinates(), 0);
241 : //cerr << "Freq " << freq << endl;
242 212 : Quantity halfpb=Quantity(1.22*C::c/freq/dishDiam_p, "rad");
243 318 : GaussianBeam newBeam(halfpb, halfpb, Quantity(0.0, "deg"));
244 107 : GaussianBeam toBeUsed;
245 : try {
246 106 : GaussianDeconvolver::deconvolve(toBeUsed, newBeam, lBeam_p);
247 : }
248 2 : catch (const AipsError& x) {
249 1 : throw(AipsError("Beam due to new effective diameter may be smaller than the beam of original dish image"));
250 : }
251 : try{
252 : //cerr <<"To be used " << toBeUsed.getMajor() << " " << toBeUsed.getMinor() <<
253 : //" " << toBeUsed.getPA() << endl;
254 105 : StokesImageUtil::Convolve(*lowIm_p, toBeUsed, true);
255 : }
256 0 : catch(const AipsError& x){
257 0 : throw(AipsError("Could not convolve SD image for some reason; try a smaller effective diameter may be"));
258 : }
259 105 : lBeam_p=newBeam;
260 : //reset cweight if it was calculated already
261 105 : cweightCalced_p=false;
262 105 : cweightApplied_p=false;
263 :
264 210 : return true;
265 : }
266 1 : void Feather::getEffectiveDishDiam(Float& xdiam, Float& ydiam){
267 1 : if(dishDiam_p <0){
268 1 : if(lBeam_p.isNull())
269 0 : throw(AipsError("No Single dish restoring beam info in image"));
270 1 : Double maj=lBeam_p.getMajor("rad");
271 1 : Double freq = worldFreq(csysHigh_p, 0);
272 : //cerr << "Freq " << freq << endl;
273 1 : dishDiam_p=1.22*C::c/freq/maj;
274 : }
275 1 : xdiam=dishDiam_p;
276 1 : ydiam=dishDiam_p;
277 1 : }
278 113 : void Feather::setSDScale(Float sdscale){
279 113 : sdScale_p=sdscale;
280 113 : }
281 :
282 0 : void Feather::getFTCutSDImage(Vector<Float>& ux, Vector<Float>& xamp, Vector<Float>& uy, Vector<Float>& yamp, const Bool radial){
283 0 : if(radial){
284 0 : uy.resize();
285 0 : yamp.resize();
286 0 : return getRadialCut(ux, xamp, *lowIm_p);
287 : }
288 :
289 0 : TempImage<Complex> cimagelow(lowIm_p->shape(), lowIm_p->coordinates() );
290 0 : StokesImageUtil::From(cimagelow, *lowIm_p);
291 0 : LatticeFFT::cfft2d( cimagelow );
292 0 : getCutXY(ux, xamp, uy, yamp, cimagelow);
293 :
294 : }
295 0 : void Feather::getFTCutIntImage(Vector<Float>& ux, Vector<Float>& xamp, Vector<Float>& uy, Vector<Float>& yamp, const Bool radial){
296 0 : if(radial){
297 0 : uy.resize();
298 0 : yamp.resize();
299 0 : return getRadialCut(ux, xamp, *highIm_p);
300 : }
301 0 : TempImage<Complex> cimagehigh(highIm_p->shape(), highIm_p->coordinates() );
302 0 : StokesImageUtil::From(cimagehigh, *highIm_p);
303 0 : LatticeFFT::cfft2d( cimagehigh );
304 0 : getCutXY(ux, xamp, uy, yamp, cimagehigh);
305 :
306 : }
307 0 : void Feather::getFeatherINT(Vector<Float>& ux, Vector<Float>& xamp, Vector<Float>& uy, Vector<Float>& yamp, Bool radial){
308 :
309 0 : calcCWeightImage();
310 0 : if(radial){
311 0 : uy.resize();
312 0 : yamp.resize();
313 0 : getRadialCut(xamp, *cwImage_p);
314 0 : getRadialUVval(xamp.shape()[0], cwImage_p->shape(), cwImage_p->coordinates(), ux);
315 :
316 :
317 : }
318 : else{
319 0 : getCutXY(ux, xamp, uy, yamp, *cwImage_p);
320 : }
321 0 : }
322 0 : void Feather::getFeatherSD(Vector<Float>& ux, Vector<Float>& xamp, Vector<Float>& uy, Vector<Float>& yamp, Bool radial, Bool normalize){
323 0 : calcCWeightImage();
324 0 : Vector<Float> xampInt, yampInt;
325 0 : if(radial){
326 0 : uy.resize();
327 0 : yamp.resize();
328 0 : getRadialCut(xampInt, *cwImage_p);
329 0 : getRadialUVval(xampInt.shape()[0], cwImage_p->shape(), cwImage_p->coordinates(), ux);
330 :
331 : }
332 : else{
333 0 : getCutXY(ux, xampInt, uy, yampInt, *cwImage_p);
334 0 : yamp.resize();
335 0 : yamp=(Float(1.0) - yampInt);
336 0 : if(!normalize)
337 0 : yamp=yamp*Float(sdScale_p*hBeam_p.getArea("arcsec2")/lBeam_p.getArea("arcsec2"));
338 : }
339 :
340 0 : xamp.resize();
341 0 : xamp=(Float(1.0) - xampInt);
342 0 : if(!normalize)
343 0 : xamp=xamp*Float(sdScale_p*hBeam_p.getArea("arcsec2")/lBeam_p.getArea("arcsec2"));
344 :
345 :
346 0 : }
347 0 : void Feather::getFeatheredCutSD(Vector<Float>& ux, Vector<Float>& xamp, Vector<Float>& uy, Vector<Float>& yamp, Bool radial){
348 :
349 0 : getFTCutSDImage(ux, xamp, uy,yamp, radial);
350 0 : xamp *=Float(sdScale_p*hBeam_p.getArea("arcsec2")/lBeam_p.getArea("arcsec2"));
351 :
352 0 : if(yamp.nelements() >0)
353 0 : yamp *=Float(sdScale_p*hBeam_p.getArea("arcsec2")/lBeam_p.getArea("arcsec2"));
354 0 : }
355 :
356 0 : void Feather::getFeatheredCutINT(Vector<Float>& ux, Vector<Float>& xamp, Vector<Float>& uy, Vector<Float>& yamp, Bool radial){
357 0 : applyFeather();
358 0 : if(radial){
359 0 : uy.resize();
360 0 : yamp.resize();
361 0 : getRadialCut(xamp, *cwHighIm_p);
362 0 : getRadialUVval(xamp.shape()[0], cwHighIm_p->shape(), cwHighIm_p->coordinates(), ux);
363 : }
364 : else
365 0 : getCutXY(ux, xamp, uy, yamp, *cwHighIm_p);
366 0 : }
367 :
368 0 : void Feather::clearWeightFlags(){
369 0 : cweightCalced_p=false;
370 0 : cweightApplied_p = false;
371 0 : }
372 :
373 112 : void Feather::applyFeather(){
374 112 : if(cweightApplied_p)
375 0 : return;
376 112 : calcCWeightImage();
377 112 : if(highIm_p.null())
378 0 : throw(AipsError("No high resolution image set"));
379 112 : cwHighIm_p=new TempImage<Complex>(highIm_p->shape(), highIm_p->coordinates() );
380 112 : StokesImageUtil::From(*cwHighIm_p, *highIm_p);
381 112 : LatticeFFT::cfft2d( *cwHighIm_p );
382 112 : Vector<Int> extraAxes(cwHighIm_p->shape().nelements()-2);
383 112 : if(extraAxes.nelements() > 0){
384 :
385 112 : if(extraAxes.nelements() ==2){
386 112 : Int n3=cwHighIm_p->shape()(2);
387 112 : Int n4=cwHighIm_p->shape()(3);
388 224 : IPosition blc(cwHighIm_p->shape());
389 224 : IPosition trc(cwHighIm_p->shape());
390 112 : blc(0)=0; blc(1)=0;
391 112 : trc(0)=cwHighIm_p->shape()(0)-1;
392 112 : trc(1)=cwHighIm_p->shape()(1)-1;
393 224 : for (Int j=0; j < n3; ++j){
394 224 : for (Int k=0; k < n4 ; ++k){
395 112 : blc(2)=j; trc(2)=j;
396 112 : blc(3)=k; trc(3)=k;
397 224 : Slicer sl(blc, trc, Slicer::endIsLast);
398 112 : SubImage<Complex> cimagehighSub(*cwHighIm_p, sl, true);
399 112 : cimagehighSub.copyData( (LatticeExpr<Complex>)((cimagehighSub * (*cwImage_p))));
400 : }
401 : }
402 : }
403 : }
404 : else{
405 0 : cwHighIm_p->copyData(
406 0 : (LatticeExpr<Complex>)(((*cwHighIm_p) * (*cwImage_p) )));
407 : }
408 112 : cweightApplied_p=true;
409 :
410 : }
411 :
412 112 : void Feather::calcCWeightImage(){
413 112 : if(cweightCalced_p)
414 0 : return;
415 112 : if(highIm_p.null())
416 0 : throw(AipsError("No Interferometer image was defined"));
417 224 : IPosition myshap(highIm_p->shape());
418 224 : Vector<Int> dirAxes;
419 112 : dirAxes=CoordinateUtil::findDirectionAxes(highIm_p->coordinates());
420 560 : for( uInt k=0; k< myshap.nelements(); ++k){
421 448 : if( (Int(k) != dirAxes[0]) && (Int(k) != dirAxes[1]))
422 224 : myshap(k)=1;
423 : }
424 :
425 112 : cwImage_p=new TempImage<Complex>(myshap, highIm_p->coordinates());
426 : {
427 224 : TempImage<Float> lowpsf(myshap, cwImage_p->coordinates());
428 112 : lowpsf.set(0.0);
429 224 : IPosition center(4, Int((myshap(0)/4)*2),
430 224 : Int((myshap(1)/4)*2),0,0);
431 112 : lowpsf.putAt(1.0, center);
432 112 : StokesImageUtil::Convolve(lowpsf, lBeam_p, false);
433 112 : StokesImageUtil::From(*cwImage_p, lowpsf);
434 : }
435 112 : LatticeFFT::cfft2d( *cwImage_p );
436 112 : LatticeExprNode node = max( *cwImage_p );
437 112 : Float fmax = abs(node.getComplex());
438 112 : cwImage_p->copyData( (LatticeExpr<Complex>)( 1.0f - (*cwImage_p)/fmax ) );
439 112 : cweightCalced_p=true;
440 112 : cweightApplied_p=false;
441 : }
442 :
443 0 : void Feather::getCutXY(Vector<Float>& ux, Vector<Float>& xamp,
444 : Vector<Float>& uy, Vector<Float>& yamp,
445 : const ImageInterface<Float>& image){
446 :
447 :
448 0 : TempImage<Complex> cimage(image.shape(), image.coordinates() );
449 0 : StokesImageUtil::From(cimage, image);
450 0 : if(image.getDefaultMask()!=""){
451 0 : ImageRegion elMask=image.getRegion(image.getDefaultMask(),
452 0 : RegionHandler::Masks);
453 0 : LatticeRegion latReg=elMask.toLatticeRegion(image.coordinates(), image.shape());
454 0 : ArrayLattice<Bool> pixmask(latReg.get());
455 0 : LatticeExpr<Complex> myexpr(iif(pixmask, cimage, Complex(0.0)) );
456 0 : cimage.copyData(myexpr);
457 :
458 : }
459 :
460 0 : LatticeFFT::cfft2d( cimage );
461 0 : getCutXY(ux, xamp, uy, yamp, cimage);
462 :
463 :
464 0 : }
465 :
466 0 : void Feather::getRadialCut(Vector<Float>& radius, Vector<Float>& radialAmp,
467 : const ImageInterface<Float>& image){
468 :
469 0 : TempImage<Complex> cimage(image.shape(), image.coordinates() );
470 0 : StokesImageUtil::From(cimage, image);
471 0 : LatticeFFT::cfft2d( cimage );
472 0 : radialAmp.resize();
473 0 : getRadialCut(radialAmp, cimage);
474 :
475 0 : getRadialUVval(radialAmp.shape()[0], cimage.shape(), cimage.coordinates(), radius);
476 : ////Get the radial value in uv- units.
477 : /* Double freq=worldFreq(image.coordinates(), 0);
478 : Int dirCoordIndex=image.coordinates().findCoordinate(Coordinate::DIRECTION);
479 : DirectionCoordinate dc=image.coordinates().directionCoordinate(dirCoordIndex);
480 : Vector<Bool> axes(2); axes(0)=true;axes(1)=true;
481 : Vector<Int> elshape(2);
482 : Vector<Int> directionIndex=CoordinateUtil::findDirectionAxes(image.coordinates());
483 : elshape(0)=image.shape()[directionIndex(0)];
484 : elshape(1)=image.shape()[directionIndex(1)];
485 : Coordinate* ftdc=dc.makeFourierCoordinate(axes,elshape);
486 : Vector<Double> xpix(radialAmp.nelements());
487 : indgen(xpix);
488 : xpix +=Double(image.shape()[directionIndex(0)])/2.0;
489 : Matrix<Double> pix(2, xpix.nelements());
490 : Matrix<Double> world(2, xpix.nelements());
491 : Vector<Bool> failures;
492 : pix.row(1)=elshape(1)/2;
493 : pix.row(0)=xpix;
494 : ftdc->toWorldMany(world, pix, failures);
495 : xpix=world.row(0);
496 : xpix=fabs(xpix)*(C::c)/freq;
497 : radius.resize(xpix.nelements());
498 : convertArray(radius, xpix);
499 : */
500 :
501 0 : }
502 0 : void Feather::getRadialUVval(const Int npix, const IPosition& imshape, const CoordinateSystem& csys, Vector<Float>& radius){
503 :
504 : ////Get the radial value in uv- units.
505 0 : Double freq=worldFreq(csys, 0);
506 0 : Int dirCoordIndex=csys.findCoordinate(Coordinate::DIRECTION);
507 0 : DirectionCoordinate dc=csys.directionCoordinate(dirCoordIndex);
508 0 : Vector<Bool> axes(2); axes(0)=true;axes(1)=true;
509 0 : Vector<Int> elshape(2);
510 0 : Vector<Int> directionIndex=CoordinateUtil::findDirectionAxes(csys);
511 0 : elshape(0)=imshape[directionIndex(0)];
512 0 : elshape(1)=imshape[directionIndex(1)];
513 0 : Coordinate* ftdc=dc.makeFourierCoordinate(axes,elshape);
514 0 : Vector<Double> xpix(npix);
515 0 : indgen(xpix);
516 0 : xpix +=Double(imshape[directionIndex(0)])/2.0;
517 0 : Matrix<Double> pix(2, npix);
518 0 : Matrix<Double> world(2, npix);
519 0 : Vector<Bool> failures;
520 0 : pix.row(1)=elshape(1)/2;
521 0 : pix.row(0)=xpix;
522 0 : ftdc->toWorldMany(world, pix, failures);
523 0 : xpix=world.row(0);
524 0 : xpix=fabs(xpix)*(C::c)/freq;
525 0 : radius.resize(xpix.nelements());
526 0 : convertArray(radius, xpix);
527 :
528 0 : }
529 0 : void Feather::getRadialCut(Vector<Float>& radialAmp, ImageInterface<Complex>& ftimage) {
530 0 : CoordinateSystem ftCoords(ftimage.coordinates());
531 : //////////////////////
532 : /*{
533 : PagedImage<Float> thisScreen(ftimage.shape(), ftCoords, "FFTimage");
534 : LatticeExpr<Float> le(abs(ftimage));
535 : thisScreen.copyData(le);
536 : }*/
537 : /////////////////////
538 0 : Vector<Int> directionIndex=CoordinateUtil::findDirectionAxes(ftCoords);
539 0 : Int spectralIndex=CoordinateUtil::findSpectralAxis(ftCoords);
540 0 : IPosition start=ftimage.shape();
541 0 : IPosition shape=ftimage.shape();
542 0 : Int centreX=shape[directionIndex(0)]/2;
543 0 : Int centreY=shape[directionIndex(1)]/2;
544 0 : Int arrLen=min(centreX, centreY);
545 :
546 0 : radialAmp.resize(arrLen);
547 0 : radialAmp.set(0.0);
548 0 : Array<Complex> tmpval;
549 :
550 0 : for(uInt k=0; k < shape.nelements(); ++k){
551 0 : start[k]=0;
552 : ///value for a single pixel in x, y, and pol
553 0 : if((k != uInt(spectralIndex)))
554 0 : shape[k]=1;
555 : }
556 : Int counter;
557 : Float sumval;
558 0 : for (Int pix =0 ; pix < arrLen; ++pix){
559 0 : sumval=0;
560 0 : counter=0;
561 0 : for (Int xval=0; xval <= pix; ++xval){
562 0 : Int yval=std::lround(sqrt(Double(pix*pix-xval*xval)));
563 0 : start[directionIndex[0]]=xval+centreX;
564 0 : start[directionIndex[1]]=yval+centreY;
565 0 : tmpval.resize();
566 0 : ftimage.getSlice(tmpval, start, shape, true);
567 0 : sumval+=fabs(mean(tmpval));
568 0 : counter+=1;
569 : }
570 0 : radialAmp[pix]=sumval/float(counter);
571 : }
572 :
573 :
574 0 : }
575 :
576 0 : void Feather::getCutXY(Vector<Float>& ux, Vector<Float>& xamp,
577 : Vector<Float>& uy, Vector<Float>& yamp, ImageInterface<Complex>& ftimage){
578 :
579 0 : CoordinateSystem ftCoords(ftimage.coordinates());
580 0 : Double freq=worldFreq(ftCoords, 0);
581 : ////
582 0 : Vector<Int> directionIndex=CoordinateUtil::findDirectionAxes(ftCoords);
583 0 : Int spectralIndex=CoordinateUtil::findSpectralAxis(ftCoords);
584 0 : Array<Complex> tmpval;
585 0 : Vector<Complex> meanval;
586 0 : IPosition start=ftimage.shape();
587 0 : IPosition shape=ftimage.shape();
588 0 : shape[directionIndex(0)]=shape[directionIndex(0)]/2;
589 0 : shape[directionIndex(1)]=shape[directionIndex(1)]/2;
590 0 : for(uInt k=0; k < shape.nelements(); ++k){
591 0 : start[k]=0;
592 0 : if((k != uInt(directionIndex(1))) && (k != uInt(spectralIndex)))
593 0 : shape[k]=1;
594 : }
595 0 : start[directionIndex(1)]=ftimage.shape()[directionIndex(1)]/2;
596 0 : start[directionIndex(0)]=ftimage.shape()[directionIndex(0)]/2;
597 0 : ftimage.getSlice(tmpval, start, shape, true);
598 0 : if(shape[spectralIndex] >1){
599 0 : meanval.resize(shape[directionIndex(1)]);
600 0 : Matrix<Complex> retmpval(tmpval);
601 0 : Bool colOrRow=spectralIndex > directionIndex(1);
602 0 : for (uInt k=0; k < meanval.nelements(); ++k){
603 0 : meanval[k]=colOrRow ? mean(retmpval.row(k)) : mean(retmpval.column(k));
604 : }
605 : }
606 : else{
607 0 : meanval=tmpval;
608 : }
609 0 : xamp.resize();
610 0 : xamp=amplitude(meanval);
611 0 : tmpval.resize();
612 0 : shape=ftimage.shape();
613 0 : shape[directionIndex(0)]=shape[directionIndex(0)]/2;
614 0 : shape[directionIndex(1)]=shape[directionIndex(1)]/2;
615 0 : for(uInt k=0; k < shape.nelements(); ++k){
616 0 : start[k]=0;
617 0 : if((k != uInt(directionIndex(0))) && (k != uInt(spectralIndex)))
618 0 : shape[k]=1;
619 : }
620 0 : start[directionIndex(1)]=ftimage.shape()[directionIndex(1)]/2;
621 0 : start[directionIndex(0)]=ftimage.shape()[directionIndex(0)]/2;
622 0 : ftimage.getSlice(tmpval, start, shape, true);
623 0 : if(shape[spectralIndex] >1){
624 0 : meanval.resize(shape[directionIndex(0)]);
625 0 : Bool colOrRow=spectralIndex > directionIndex(0);
626 0 : Matrix<Complex> retmpval(tmpval);
627 0 : for (uInt k=0; k < meanval.nelements(); ++k){
628 0 : meanval[k]=colOrRow ? mean(retmpval.row(k)) : mean(retmpval.column(k));
629 : }
630 : }
631 : else{
632 0 : meanval.resize( tmpval.size());
633 0 : meanval=tmpval;
634 : }
635 0 : yamp.resize();
636 0 : yamp=amplitude(meanval);
637 0 : Int dirCoordIndex=ftCoords.findCoordinate(Coordinate::DIRECTION);
638 0 : DirectionCoordinate dc=ftCoords.directionCoordinate(dirCoordIndex);
639 0 : Vector<Bool> axes(2); axes(0)=true;axes(1)=true;
640 0 : Vector<Int> elshape(2);
641 0 : elshape(0)=ftimage.shape()[directionIndex(0)];
642 0 : elshape(1)=ftimage.shape()[directionIndex(1)];
643 0 : Coordinate* ftdc=dc.makeFourierCoordinate(axes,elshape);
644 0 : Vector<Double> xpix(xamp.nelements());
645 0 : indgen(xpix);
646 0 : xpix +=Double(ftimage.shape()[directionIndex(0)])/2.0;
647 0 : Matrix<Double> pix(2, xpix.nelements());
648 0 : Matrix<Double> world(2, xpix.nelements());
649 0 : Vector<Bool> failures;
650 0 : pix.row(1)=elshape(0)/2;
651 0 : pix.row(0)=xpix;
652 0 : ftdc->toWorldMany(world, pix, failures);
653 0 : xpix=world.row(0);
654 : //cerr << "xpix " << xpix << endl;
655 0 : xpix=fabs(xpix)*(C::c)/freq;
656 0 : Vector<Double> ypix(yamp.nelements());
657 0 : indgen(ypix);
658 0 : ypix +=Double(ftimage.shape()[directionIndex(1)])/2.0;
659 0 : pix.resize(2, ypix.nelements());
660 0 : world.resize();
661 0 : pix.row(1)=ypix;
662 0 : pix.row(0)=elshape(1)/2;
663 0 : ftdc->toWorldMany(world, pix, failures);
664 0 : ypix=world.row(1);
665 0 : ypix=fabs(ypix)*(C::c/freq);
666 0 : ux.resize(xpix.nelements());
667 0 : uy.resize(ypix.nelements());
668 0 : convertArray(ux, xpix);
669 0 : convertArray(uy, ypix);
670 :
671 :
672 0 : }
673 :
674 :
675 :
676 112 : Bool Feather::saveFeatheredImage(const String& imagename){
677 112 : applyFeather();
678 : Int stokesAxis, spectralAxis;
679 112 : spectralAxis=CoordinateUtil::findSpectralAxis(csysHigh_p);
680 224 : Vector<Stokes::StokesTypes> stokesvec;
681 112 : stokesAxis=CoordinateUtil::findStokesAxis(stokesvec, csysHigh_p);
682 224 : Vector<Int> dirAxes=CoordinateUtil::findDirectionAxes(csysHigh_p);
683 112 : Int n3=highIm_p->shape()(stokesAxis);
684 112 : Int n4=highIm_p->shape()(spectralAxis);
685 224 : IPosition blc(highIm_p->shape());
686 224 : IPosition trc(highIm_p->shape());
687 112 : blc(dirAxes(0))=0; blc(dirAxes(1))=0;
688 112 : trc(dirAxes(0))=highIm_p->shape()(dirAxes(0))-1;
689 112 : trc(dirAxes(1))=highIm_p->shape()(dirAxes(1))-1;
690 336 : TempImage<Complex>cimagelow(lowIm_p->shape(), lowIm_p->coordinates());
691 112 : StokesImageUtil::From(cimagelow, *lowIm_p);
692 112 : LatticeFFT::cfft2d( cimagelow);
693 112 : Float sdScaling = sdScale_p*hBeam_p.getArea("arcsec2")/lBeam_p.getArea("arcsec2");
694 224 : for (Int j=0; j < n3; ++j){
695 224 : for (Int k=0; k < n4 ; ++k){
696 112 : blc(stokesAxis)=j; trc(stokesAxis)=j;
697 112 : blc(spectralAxis)=k; trc(spectralAxis)=k;
698 224 : Slicer sl(blc, trc, Slicer::endIsLast);
699 224 : SubImage<Complex> cimagehighSub(*cwHighIm_p, sl, true);
700 112 : SubImage<Complex> cimagelowSub(cimagelow, sl, true);
701 112 : cimagelowSub.copyData( (LatticeExpr<Complex>)((cimagehighSub + cimagelowSub * sdScaling)));
702 : }
703 : }
704 : // FT back to image plane
705 112 : LatticeFFT::cfft2d( cimagelow, false);
706 :
707 : // write to output image
708 336 : PagedImage<Float> featherImage(highIm_p->shape(), highIm_p->coordinates(), imagename );
709 112 : StokesImageUtil::To(featherImage, cimagelow);
710 112 : ImageUtilities::copyMiscellaneous(featherImage, *highIm_p);
711 224 : String maskofHigh=highIm_p->getDefaultMask();
712 112 : String maskofLow=lowIm_p->getDefaultMask();
713 112 : if(maskofHigh != String("")){
714 : //ImageUtilities::copyMask(featherImage, *highIm_p, "mask0", maskofHigh, AxesSpecifier());
715 : //featherImage.setDefaultMask("mask0");
716 16 : Imager::copyMask(featherImage, *highIm_p, maskofHigh);
717 16 : if(maskofLow != String("")){
718 16 : featherImage.pixelMask().copyData(LatticeExpr<Bool> (featherImage.pixelMask() && (lowIm_p->pixelMask())));
719 :
720 : }
721 : }
722 96 : else if(maskofLow != String("")){
723 : //ImageUtilities::copyMask(featherImage, *lowIm_p, "mask0", maskofLow, AxesSpecifier());
724 16 : Imager::copyMask(featherImage, *lowIm_p, maskofLow);
725 : }
726 :
727 224 : return true;
728 : }
729 :
730 0 : void Feather::applyDishDiam(ImageInterface<Complex>& image, GaussianBeam& beam, Float effDiam, ImageInterface<Float>& fftim, Vector<Quantity>& extraconv){
731 : /*
732 : MathFunc<Float> dd(SPHEROIDAL);
733 : Vector<Float> valsph(31);
734 : for(Int k=0; k <31; ++k){
735 : valsph(k)=dd.value((Float)(k)/10.0);
736 : }
737 : Quantity fulrad((52.3/2.0*3.0/1.1117),"arcsec");
738 : Quantity lefreq(224.0/effDiam*9.0, "GHz");
739 :
740 : PBMath1DNumeric elpb(valsph, fulrad, lefreq);*/
741 : ////////////////////////
742 : /*{
743 : PagedImage<Float> thisScreen(image.shape(), image.coordinates(), "Before_apply.image");
744 : LatticeExpr<Float> le(abs(image));
745 : thisScreen.copyData(le);
746 : }*/
747 : //////////////////////
748 :
749 0 : Double freq = worldFreq(image.coordinates(), 0);
750 : //cerr << "Freq " << freq << endl;
751 0 : Quantity halfpb=Quantity(1.22*C::c/freq/effDiam, "rad");
752 : //PBMath1DAiry elpb( Quantity(effDiam,"m"), Quantity(0.01,"m"),
753 : // Quantity(0.8564,"deg"), Quantity(1.0,"GHz"));
754 0 : GaussianBeam newBeam(halfpb, halfpb, Quantity(0.0, "deg"));
755 : try {
756 :
757 0 : GaussianBeam toBeUsed;
758 : //cerr << "beam " << beam.toVector() << endl;
759 : //cerr << "newBeam " << newBeam.toVector() << endl;
760 0 : GaussianDeconvolver::deconvolve(toBeUsed, newBeam, beam);
761 0 : extraconv.resize(3);
762 : // use the Major difference
763 0 : extraconv(0) = toBeUsed.getMajor();
764 0 : extraconv(1) = toBeUsed.getMinor();
765 0 : extraconv(2) = toBeUsed.getPA();
766 :
767 : }
768 0 : catch (const AipsError& x) {
769 0 : throw(AipsError("Beam due to new effective diameter may be smaller than the beam of original dish image"));
770 : }
771 : //////////////////////
772 : //1 GHz equivalent
773 : // halfpb=Quantity(1.22*C::c/1.0e9/effDiam, "rad");
774 : //cerr << "halfpb " << halfpb << endl;
775 : //PBMath1DGauss elpb(halfpb, Quantity(0.8564,"deg"), Quantity(1.0,"GHz"), true);
776 :
777 0 : fftim.set(0.0);
778 :
779 0 : IPosition center(4, Int((fftim.shape()(0)/4)*2),
780 0 : Int((fftim.shape()(1)/4)*2),0,0);
781 0 : fftim.putAt(1.0, center);
782 0 : StokesImageUtil::Convolve(fftim, newBeam, false);
783 0 : StokesImageUtil::From(image, fftim);
784 : /*
785 : TempImage<Complex> elbeamo(image.shape(), image.coordinates());
786 : elbeamo.set(1.0);
787 : MDirection wcenter;
788 : {
789 : Int directionIndex=
790 : image.coordinates().findCoordinate(Coordinate::DIRECTION);
791 : DirectionCoordinate
792 : directionCoord=image.coordinates().directionCoordinate(directionIndex);
793 : Vector<Double> pcenter(2);
794 : pcenter(0) = image.shape()(directionIndex)/2;
795 : pcenter(1) = image.shape()(directionIndex+1)/2;
796 : directionCoord.toWorld( wcenter, pcenter );
797 : }
798 : elpb.applyPB(elbeamo, elbeamo, wcenter, Quantity(0.0, "deg"),
799 : BeamSquint::NONE);
800 :
801 :
802 : StokesImageUtil::To(fftim, elbeamo);
803 : */
804 : ///////////////////
805 : //StokesImageUtil::FitGaussianPSF(fftim,
806 : // beam);
807 : //cerr << "To be convolved beam 2 " << beam.toVector() << endl;
808 : ////////////////
809 :
810 0 : LatticeFFT::cfft2d(image);
811 : //image.copyData((LatticeExpr<Complex>)(elbeamo) );
812 : //elbeamo.copyData(image);
813 : // LatticeFFT::cfft2d(elbeamo, false);
814 : //StokesImageUtil::To(fftim, elbeamo);
815 : //StokesImageUtil::FitGaussianPSF(fftim,
816 : // beam);
817 : /////////
818 : /*{
819 : PagedImage<Float> thisScreen(image.shape(), image.coordinates(), "After_apply.image");
820 : //LatticeExpr<Float> le(abs(image));
821 : thisScreen.copyData(fftim);
822 : }*/
823 : ////////
824 0 : }
825 :
826 :
827 : /*
828 : void Feather::feather(const String& image, const ImageInterface<Float>& high, const ImageInterface<Float>& low0, const Float& sdScale, const String& lowPSF, const Bool useDefaultPB, const String& vpTableStr, Float effDiam, const Bool doPlot){
829 :
830 : LogIO os(LogOrigin("Feather", "feather()", WHERE));
831 : try {
832 :
833 :
834 : GaussianBeam hBeam , lBeam;
835 : Vector<Quantity> extraconv;
836 : ImageInfo highInfo=high.imageInfo();
837 : hBeam=highInfo.restoringBeam();
838 : ImageInfo lowInfo=low0.imageInfo();
839 : lBeam=lowInfo.restoringBeam();
840 : if(hBeam.isNull()) {
841 : os << LogIO::WARN
842 : << "High resolution image does not have any resolution information - will be unable to scale correctly.\n"
843 : << LogIO::POST;
844 : }
845 :
846 : PBMath * myPBp = 0;
847 : if((lowPSF=="") && lBeam.isNull()) {
848 : // create the low res's PBMath object, needed to apply PB
849 : // to make high res Fourier weight image
850 : if (useDefaultPB) {
851 : // look up the telescope in ObsInfo
852 : ObsInfo oi = low0.coordinates().obsInfo();
853 : String myTelescope = oi.telescope();
854 : if (myTelescope == "") {
855 : os << LogIO::SEVERE << "No telescope imbedded in low res image"
856 : << LogIO::POST;
857 : os << LogIO::SEVERE << "Create a PB description with the vpmanager"
858 : << LogIO::EXCEPTION;
859 : }
860 : Quantity qFreq;
861 : {
862 : Double freq = worldFreq(low0.coordinates(), 0);
863 : qFreq = Quantity( freq, "Hz" );
864 : }
865 : String band;
866 : PBMath::CommonPB whichPB;
867 : String pbName;
868 : // get freq from coordinates
869 : PBMath::whichCommonPBtoUse (myTelescope, qFreq, band, whichPB,
870 : pbName);
871 : if (whichPB == PBMath::UNKNOWN) {
872 : os << LogIO::SEVERE << "Unknown telescope for PB type: "
873 : << myTelescope << LogIO::EXCEPTION;
874 : }
875 : myPBp = new PBMath(whichPB);
876 : } else {
877 : // get the PB from the vpTable
878 : Table vpTable( vpTableStr );
879 : ScalarColumn<TableRecord> recCol(vpTable, (String)"pbdescription");
880 : myPBp = new PBMath(recCol(0));
881 : }
882 : AlwaysAssert((myPBp != 0), AipsError);
883 : }
884 :
885 : // regrid the single dish image
886 : TempImage<Float> low(high.shape(), high.coordinates());
887 : {
888 : IPosition axes(2,0,1);
889 : if(high.shape().nelements() >2){
890 : Int spectralAxisIndex=high.coordinates().
891 : findCoordinate(Coordinate::SPECTRAL);
892 : if(spectralAxisIndex > -1){
893 : axes.resize(3);
894 : axes(0)=0;
895 : axes(1)=1;
896 : axes(2)=spectralAxisIndex+1;
897 : }
898 : }
899 : ImageRegrid<Float> ir;
900 : ir.regrid(low, Interpolate2D::LINEAR, axes, low0);
901 : }
902 :
903 : // get image center direction (needed for SD PB, which is needed for
904 : // the high res Fourier weight image
905 : MDirection wcenter;
906 : {
907 : Int directionIndex=
908 : high.coordinates().findCoordinate(Coordinate::DIRECTION);
909 : AlwaysAssert(directionIndex>=0, AipsError);
910 : DirectionCoordinate
911 : directionCoord=high.coordinates().directionCoordinate(directionIndex);
912 : Vector<Double> pcenter(2);
913 : pcenter(0) = high.shape()(0)/2;
914 : pcenter(1) = high.shape()(1)/2;
915 : directionCoord.toWorld( wcenter, pcenter );
916 : }
917 :
918 : // make the weight image for high res Fourier plane: 1 - normalized(FT(sd_PB))
919 : IPosition myshap(high.shape());
920 : for( uInt k=2; k< myshap.nelements(); ++k){
921 : myshap(k)=1;
922 : }
923 :
924 : TempImage<Float> lowpsf0;
925 : TempImage<Complex> cweight(myshap, high.coordinates());
926 : if(lowPSF=="") {
927 : os << LogIO::NORMAL // Loglevel INFO
928 : << "Using primary beam to determine weighting.\n" << LogIO::POST;
929 : if(lBeam.isNull()) {
930 : cweight.set(1.0);
931 : if (myPBp != 0) {
932 : myPBp->applyPB(cweight, cweight, wcenter, Quantity(0.0, "deg"),
933 : BeamSquint::NONE);
934 :
935 : lowpsf0=TempImage<Float>(cweight.shape(), cweight.coordinates());
936 :
937 : os << LogIO::NORMAL // Loglevel INFO
938 : << "Determining scaling from SD Primary Beam.\n"
939 : << LogIO::POST;
940 : StokesImageUtil::To(lowpsf0, cweight);
941 : StokesImageUtil::FitGaussianPSF(lowpsf0,
942 : lBeam);
943 : }
944 : delete myPBp;
945 : }
946 : else{
947 : os << LogIO::NORMAL // Loglevel INFO
948 : << "Determining scaling from SD restoring beam.\n"
949 : << LogIO::POST;
950 : lowpsf0=TempImage<Float>(cweight.shape(), cweight.coordinates());
951 : lowpsf0.set(0.0);
952 : IPosition center(4, Int((cweight.shape()(0)/4)*2),
953 : Int((cweight.shape()(1)/4)*2),0,0);
954 : lowpsf0.putAt(1.0, center);
955 : StokesImageUtil::Convolve(lowpsf0, lBeam, false);
956 : StokesImageUtil::From(cweight, lowpsf0);
957 :
958 : }
959 : }
960 : else {
961 : os << LogIO::NORMAL // Loglevel INFO
962 : << "Using specified low resolution PSF to determine weighting.\n"
963 : << LogIO::POST;
964 : // regrid the single dish psf
965 : PagedImage<Float> lowpsfDisk(lowPSF);
966 : IPosition lshape(lowpsfDisk.shape());
967 : lshape.resize(4);
968 : lshape(2)=1; lshape(3)=1;
969 : TempImage<Float>lowpsf(lshape,lowpsfDisk.coordinates());
970 : IPosition blc(lowpsfDisk.shape());
971 : IPosition trc(lowpsfDisk.shape());
972 : blc(0)=0; blc(1)=0;
973 : trc(0)=lowpsfDisk.shape()(0)-1;
974 : trc(1)=lowpsfDisk.shape()(1)-1;
975 : for( uInt k=2; k < lowpsfDisk.shape().nelements(); ++k){
976 : blc(k)=0; trc(k)=0;
977 : }// taking first plane
978 : Slicer sl(blc, trc, Slicer::endIsLast);
979 : lowpsf.copyData(SubImage<Float>(lowpsfDisk, sl, false));
980 : lowpsf0=TempImage<Float> (myshap, high.coordinates());
981 : {
982 : ImageRegrid<Float> ir;
983 : IPosition axes(2,0,1); // if its a cube, regrid the spectral too
984 : ir.regrid(lowpsf0, Interpolate2D::LINEAR, axes, lowpsf);
985 : }
986 : if(lBeam.isNull()) {
987 : os << LogIO::NORMAL // Loglevel INFO
988 : << "Determining scaling from low resolution PSF.\n" << LogIO::POST;
989 : StokesImageUtil::FitGaussianPSF(lowpsf, lBeam);
990 : }
991 : StokesImageUtil::From(cweight, lowpsf0);
992 : }
993 :
994 : LatticeFFT::cfft2d( cweight );
995 : if(effDiam >0.0){
996 : //cerr << "in effdiam" << effDiam << endl;
997 : applyDishDiam(cweight, lBeam, effDiam, lowpsf0, extraconv);
998 : lowpsf0.copyData((LatticeExpr<Float>)(lowpsf0/max(lowpsf0)));
999 : StokesImageUtil::FitGaussianPSF(lowpsf0, lBeam);
1000 :
1001 : Int directionIndex=
1002 : cweight.coordinates().findCoordinate(Coordinate::DIRECTION);
1003 : Image2DConvolver<Float>::convolve(
1004 : os, low, low, VectorKernel::toKernelType("gauss"), IPosition(2, directionIndex, directionIndex+1),
1005 : extraconv, true, 1.0, true
1006 : );
1007 :
1008 : }
1009 : LatticeExprNode node = max( cweight );
1010 : Float fmax = abs(node.getComplex());
1011 : cweight.copyData( (LatticeExpr<Complex>)( 1.0f - cweight/fmax ) );
1012 : //Plotting part
1013 : if(doPlot){
1014 : CoordinateSystem ftCoords(cweight.coordinates());
1015 : Double freq=worldFreq(ftCoords, 0);
1016 : ////
1017 : Int directionIndex=ftCoords.findCoordinate(Coordinate::DIRECTION);
1018 : Array<Complex> tmpval;
1019 : IPosition start=cweight.shape();
1020 : IPosition shape=cweight.shape()/2;
1021 : for(uInt k=0; k < shape.nelements(); ++k){
1022 : start[k]=0;
1023 : if(k != uInt(directionIndex+1))
1024 : shape[k]=1;
1025 : }
1026 : start[directionIndex+1]=cweight.shape()[directionIndex+1]/2;
1027 : start[directionIndex]=cweight.shape()[directionIndex]/2;
1028 : cweight.getSlice(tmpval, start, shape, true);
1029 : Vector<Float> x=amplitude(tmpval);
1030 : Vector<Float> xdish=(Float(1.0) - x)*Float(sdScale);
1031 : tmpval.resize();
1032 : shape=cweight.shape()/2;
1033 : for(uInt k=0; k < shape.nelements(); ++k){
1034 : start[k]=0;
1035 : if(k != uInt(directionIndex))
1036 : shape[k]=1;
1037 : }
1038 : start[directionIndex+1]=cweight.shape()[directionIndex+1]/2;
1039 : start[directionIndex]=cweight.shape()[directionIndex]/2;
1040 : cweight.getSlice(tmpval, start, shape, true);
1041 : Vector<Float> y=amplitude(tmpval);
1042 : Vector<Float> ydish=(Float(1.0)-y)*Float(sdScale);
1043 : DirectionCoordinate dc=ftCoords.directionCoordinate(directionIndex);
1044 : Vector<Bool> axes(2); axes(0)=true;axes(1)=true;
1045 : Vector<Int> elshape(2);
1046 : elshape(0)=cweight.shape()[directionIndex];
1047 : elshape(1)=cweight.shape()[directionIndex+1];
1048 : Coordinate* ftdc=dc.makeFourierCoordinate(axes,elshape);
1049 : Vector<Double> xpix(x.nelements());
1050 : indgen(xpix);
1051 : xpix +=Double(cweight.shape()[directionIndex])/2.0;
1052 : Matrix<Double> pix(2, xpix.nelements());
1053 : Matrix<Double> world(2, xpix.nelements());
1054 : Vector<Bool> failures;
1055 : pix.row(1)=elshape(0)/2;
1056 : pix.row(0)=xpix;
1057 : ftdc->toWorldMany(world, pix, failures);
1058 : xpix=world.row(0);
1059 : //cerr << "xpix " << xpix << endl;
1060 : xpix=fabs(xpix)*(C::c)/freq;
1061 : Vector<Double> ypix(y.nelements());
1062 : indgen(ypix);
1063 : ypix +=Double(cweight.shape()[directionIndex+1])/2.0;
1064 : pix.resize(2, ypix.nelements());
1065 : world.resize();
1066 : pix.row(1)=ypix;
1067 : pix.row(0)=elshape(1)/2;
1068 : ftdc->toWorldMany(world, pix, failures);
1069 : ypix=world.row(1);
1070 : ypix=fabs(ypix)*(C::c/freq);
1071 : PlotServerProxy *plotter = dbus::launch<PlotServerProxy>( );
1072 : dbus::variant panel_id = plotter->panel( "Feather function slice cuts", "Distance in m", "Amp", "Feather",
1073 : std::vector<int>( ), "right");
1074 :
1075 : if ( panel_id.type( ) != dbus::variant::INT ) {
1076 : os << LogIO::SEVERE << "failed to start plotter" << LogIO::POST;
1077 : return ;
1078 : }
1079 :
1080 :
1081 : plotter->line(dbus::af(xpix),dbus::af(x),"blue","x-interferometer weight",panel_id.getInt( ));
1082 : plotter->line(dbus::af(xpix),dbus::af(xdish),"cyan","x-singledish weight",panel_id.getInt( ));
1083 :
1084 : plotter->line(dbus::af(ypix),dbus::af(y),"green","y-interferometer weight",panel_id.getInt( ));
1085 : plotter->line(dbus::af(ypix),dbus::af(ydish),"purple","y-singledish weight",panel_id.getInt( ));
1086 :
1087 : }
1088 : //End plotting part
1089 : // FT high res image
1090 : TempImage<Complex> cimagehigh(high.shape(), high.coordinates() );
1091 : StokesImageUtil::From(cimagehigh, high);
1092 : LatticeFFT::cfft2d( cimagehigh );
1093 :
1094 : // FT low res image
1095 : TempImage<Complex> cimagelow(high.shape(), high.coordinates() );
1096 : StokesImageUtil::From(cimagelow, low);
1097 : LatticeFFT::cfft2d( cimagelow );
1098 :
1099 :
1100 : // This factor comes from the beam volumes
1101 : if(sdScale != 1.0)
1102 : os << LogIO::NORMAL // Loglevel INFO
1103 : << "Multiplying single dish data by user specified factor "
1104 : << sdScale << ".\n" << LogIO::POST;
1105 : Float sdScaling = sdScale;
1106 : if (! hBeam.isNull() && ! lBeam.isNull()) {
1107 :
1108 : Float beamFactor=
1109 : hBeam.getArea("arcsec2")/lBeam.getArea("arcsec2");
1110 : os << LogIO::NORMAL // Loglevel INFO
1111 : << "Applying additional scaling for ratio of the volumes of the high to the low resolution images : "
1112 : << beamFactor << ".\n" << LogIO::POST;
1113 : sdScaling*=beamFactor;
1114 : }
1115 : else {
1116 : os << LogIO::WARN << "Insufficient information to scale correctly.\n"
1117 : << LogIO::POST;
1118 : }
1119 :
1120 : // combine high and low res, appropriately normalized, in Fourier
1121 : // plane. The vital point to remember is that cimagelow is already
1122 : // multiplied by 1-cweight so we only need adjust for the ratio of beam
1123 : // volumes
1124 : Vector<Int> extraAxes(cimagehigh.shape().nelements()-2);
1125 : if(extraAxes.nelements() > 0){
1126 :
1127 : if(extraAxes.nelements() ==2){
1128 : Int n3=cimagehigh.shape()(2);
1129 : Int n4=cimagehigh.shape()(3);
1130 : IPosition blc(cimagehigh.shape());
1131 : IPosition trc(cimagehigh.shape());
1132 : blc(0)=0; blc(1)=0;
1133 : trc(0)=cimagehigh.shape()(0)-1;
1134 : trc(1)=cimagehigh.shape()(1)-1;
1135 : for (Int j=0; j < n3; ++j){
1136 : for (Int k=0; k < n4 ; ++k){
1137 : blc(2)=j; trc(2)=j;
1138 : blc(3)=k; trc(3)=k;
1139 : Slicer sl(blc, trc, Slicer::endIsLast);
1140 : SubImage<Complex> cimagehighSub(cimagehigh, sl, true);
1141 : SubImage<Complex> cimagelowSub(cimagelow, sl, true);
1142 : cimagehighSub.copyData( (LatticeExpr<Complex>)((cimagehighSub * cweight + cimagelowSub * sdScaling)));
1143 : }
1144 : }
1145 : }
1146 : }
1147 : else{
1148 : cimagehigh.copyData(
1149 : (LatticeExpr<Complex>)((cimagehigh * cweight
1150 : + cimagelow * sdScaling)));
1151 : }
1152 : // FT back to image plane
1153 : LatticeFFT::cfft2d( cimagehigh, false);
1154 :
1155 : // write to output image
1156 : PagedImage<Float> featherImage(high.shape(), high.coordinates(), image );
1157 : StokesImageUtil::To(featherImage, cimagehigh);
1158 : ImageUtilities::copyMiscellaneous(featherImage, high);
1159 :
1160 : try{
1161 : { // write data processing history into image logtable
1162 : LoggerHolder imagelog (false);
1163 : LogSink& sink = imagelog.sink();
1164 : LogOrigin lor(String("Feather"), String("feather()"));
1165 : LogMessage msg(lor);
1166 : ostringstream oos;
1167 : oos << endl << "Imager::feather() input paramaters:" << endl
1168 : << "Feathered image = '" << image << "'" << endl
1169 : << "High resolution image ='" << high.name() << "'" << endl
1170 : << "Low resolution image = '" << low0.name() << "'" << endl
1171 : << "Low resolution PSF = '" << lowPSF << "'" << endl << endl;
1172 : String inputs(oos);
1173 : sink.postLocally(msg.message(inputs));
1174 : imagelog.flush();
1175 :
1176 : LoggerHolder& log = featherImage.logger();
1177 : log.append(imagelog);
1178 : log.flush();
1179 : }
1180 : }catch(exception& x){
1181 :
1182 : os << LogIO::WARN << "Caught exception: " << x.what()
1183 : << LogIO::POST;
1184 :
1185 :
1186 : }
1187 : catch(...){
1188 : os << LogIO::SEVERE << "Unknown exception handled" << LogIO::POST;
1189 :
1190 : }
1191 :
1192 :
1193 :
1194 :
1195 :
1196 :
1197 :
1198 :
1199 : }catch(exception& x){
1200 :
1201 : os << LogIO::WARN << "Caught exception: " << x.what()
1202 : << LogIO::POST;
1203 : }
1204 :
1205 :
1206 :
1207 : }
1208 : */
1209 113 : void Feather::feather(const String& image, const ImageInterface<Float>& high, const ImageInterface<Float>& low0, const Float& sdScale, const String& lowPSF, const Bool useDefaultPB, const String& vpTableStr, Float effDiam, const Bool doHiPassFilterOnSD){
1210 :
1211 226 : Feather plume;
1212 113 : plume.setINTImage(high);
1213 226 : ImageInfo lowInfo=low0.imageInfo();
1214 226 : GaussianBeam lBeam=lowInfo.restoringBeam();
1215 113 : if(lBeam.isNull()) {
1216 0 : getLowBeam(low0, lowPSF, useDefaultPB, vpTableStr, lBeam);
1217 0 : if(lBeam.isNull())
1218 0 : throw(AipsError("Do not have low resolution beam info "));
1219 0 : TempImage<Float> newLow(low0.shape(), low0.coordinates());
1220 0 : newLow.copyData(low0);
1221 0 : lowInfo.removeRestoringBeam();
1222 0 : lowInfo.setRestoringBeam(lBeam);
1223 0 : newLow.setImageInfo(lowInfo);
1224 0 : plume.setSDImage(newLow);
1225 : }
1226 : else{
1227 113 : plume.setSDImage(low0);
1228 : }
1229 113 : plume.setSDScale(sdScale);
1230 113 : if(effDiam >0.0){
1231 105 : plume.setEffectiveDishDiam(effDiam, effDiam);
1232 : }
1233 8 : else if(doHiPassFilterOnSD){
1234 : Float xdiam, ydiam;
1235 1 : plume.getEffectiveDishDiam(xdiam, ydiam);
1236 1 : plume.setEffectiveDishDiam(xdiam, ydiam);
1237 : }
1238 112 : plume.saveFeatheredImage(image);
1239 :
1240 112 : }
1241 :
1242 :
1243 0 : void Feather::getLowBeam(const ImageInterface<Float>& low0, const String& lowPSF, const Bool useDefaultPB, const String& vpTableStr, GaussianBeam& lBeam){
1244 0 : LogIO os(LogOrigin("Feather", "getLowBeam()", WHERE));
1245 0 : PBMath * myPBp = 0;
1246 0 : if((lowPSF=="") && lBeam.isNull()) {
1247 : // create the low res's PBMath object, needed to apply PB
1248 : // to make high res Fourier weight image
1249 0 : if (useDefaultPB) {
1250 : // look up the telescope in ObsInfo
1251 0 : ObsInfo oi = low0.coordinates().obsInfo();
1252 0 : String myTelescope = oi.telescope();
1253 0 : if (myTelescope == "") {
1254 : os << LogIO::SEVERE << "No telescope imbedded in low res image"
1255 0 : << LogIO::POST;
1256 : os << LogIO::SEVERE << "Create a PB description with the vpmanager"
1257 0 : << LogIO::EXCEPTION;
1258 : }
1259 0 : Quantity qFreq;
1260 : {
1261 0 : Double freq = worldFreq(low0.coordinates(), 0);
1262 0 : qFreq = Quantity( freq, "Hz" );
1263 : }
1264 0 : String band;
1265 : PBMath::CommonPB whichPB;
1266 0 : String pbName;
1267 : // get freq from coordinates
1268 0 : PBMath::whichCommonPBtoUse (myTelescope, qFreq, band, whichPB,
1269 : pbName);
1270 0 : if (whichPB == PBMath::UNKNOWN) {
1271 : os << LogIO::SEVERE << "Unknown telescope for PB type: "
1272 0 : << myTelescope << LogIO::EXCEPTION;
1273 : }
1274 0 : myPBp = new PBMath(whichPB);
1275 : } else {
1276 : // get the PB from the vpTable
1277 0 : Table vpTable( vpTableStr );
1278 0 : ScalarColumn<TableRecord> recCol(vpTable, (String)"pbdescription");
1279 0 : myPBp = new PBMath(recCol(0));
1280 : }
1281 0 : AlwaysAssert((myPBp != 0), AipsError);
1282 : }
1283 :
1284 :
1285 : // get image center direction
1286 0 : MDirection wcenter;
1287 : {
1288 : Int directionIndex=
1289 0 : low0.coordinates().findCoordinate(Coordinate::DIRECTION);
1290 0 : AlwaysAssert(directionIndex>=0, AipsError);
1291 : DirectionCoordinate
1292 0 : directionCoord=low0.coordinates().directionCoordinate(directionIndex);
1293 0 : Vector<Double> pcenter(2);
1294 0 : pcenter(0) = low0.shape()(0)/2;
1295 0 : pcenter(1) = low0.shape()(1)/2;
1296 0 : directionCoord.toWorld( wcenter, pcenter );
1297 : }
1298 :
1299 : // make the weight image
1300 0 : IPosition myshap(low0.shape());
1301 0 : for( uInt k=2; k< myshap.nelements(); ++k){
1302 0 : myshap(k)=1;
1303 : }
1304 :
1305 0 : TempImage<Float> lowpsf0;
1306 0 : TempImage<Complex> cweight(myshap, low0.coordinates());
1307 0 : if(lowPSF=="") {
1308 : os << LogIO::NORMAL // Loglevel INFO
1309 0 : << "Using primary beam to determine weighting.\n" << LogIO::POST;
1310 0 : cweight.set(1.0);
1311 0 : if (myPBp != 0) {
1312 0 : myPBp->applyPB(cweight, cweight, wcenter, Quantity(0.0, "deg"),
1313 0 : BeamSquint::NONE);
1314 :
1315 0 : lowpsf0=TempImage<Float>(cweight.shape(), cweight.coordinates());
1316 :
1317 : os << LogIO::NORMAL // Loglevel INFO
1318 : << "Determining scaling from SD Primary Beam.\n"
1319 0 : << LogIO::POST;
1320 0 : StokesImageUtil::To(lowpsf0, cweight);
1321 0 : StokesImageUtil::FitGaussianPSF(lowpsf0,
1322 : lBeam);
1323 : }
1324 0 : delete myPBp;
1325 : }
1326 : else {
1327 : os << LogIO::NORMAL // Loglevel INFO
1328 : << "Using specified low resolution PSF to determine weighting.\n"
1329 0 : << LogIO::POST;
1330 : // regrid the single dish psf
1331 0 : PagedImage<Float> lowpsfDisk(lowPSF);
1332 0 : IPosition lshape(lowpsfDisk.shape());
1333 0 : lshape.resize(4);
1334 0 : lshape(2)=1; lshape(3)=1;
1335 0 : TempImage<Float>lowpsf(lshape,lowpsfDisk.coordinates());
1336 0 : IPosition blc(lowpsfDisk.shape());
1337 0 : IPosition trc(lowpsfDisk.shape());
1338 0 : blc(0)=0; blc(1)=0;
1339 0 : trc(0)=lowpsfDisk.shape()(0)-1;
1340 0 : trc(1)=lowpsfDisk.shape()(1)-1;
1341 0 : for( uInt k=2; k < lowpsfDisk.shape().nelements(); ++k){
1342 0 : blc(k)=0; trc(k)=0;
1343 : }// taking first plane
1344 0 : Slicer sl(blc, trc, Slicer::endIsLast);
1345 0 : lowpsf.copyData(SubImage<Float>(lowpsfDisk, sl, false));
1346 : os << LogIO::NORMAL // Loglevel INFO
1347 0 : << "Determining scaling from low resolution PSF.\n" << LogIO::POST;
1348 0 : StokesImageUtil::FitGaussianPSF(lowpsf, lBeam);
1349 : }
1350 :
1351 :
1352 :
1353 0 : }
1354 :
1355 107 : Double Feather::worldFreq(const CoordinateSystem& cs, Int spectralpix){
1356 : ///freq part
1357 107 : Int spectralIndex=cs.findCoordinate(Coordinate::SPECTRAL);
1358 : SpectralCoordinate
1359 : spectralCoord=
1360 214 : cs.spectralCoordinate(spectralIndex);
1361 214 : Vector<String> units(1); units = "Hz";
1362 107 : spectralCoord.setWorldAxisUnits(units);
1363 214 : Vector<Double> spectralWorld(1);
1364 107 : Vector<Double> spectralPixel(1);
1365 107 : spectralPixel(0) = spectralpix;
1366 107 : spectralCoord.toWorld(spectralWorld, spectralPixel);
1367 107 : Double freq = spectralWorld(0);
1368 214 : return freq;
1369 : }
1370 :
1371 :
1372 :
1373 :
1374 :
1375 : }//# NAMESPACE CASA - END
|
