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