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 17 : Feather::Feather(): dishDiam_p(-1.0), cweightCalced_p(false), cweightApplied_p(false), sdScale_p(1.0){
71 17 : highIm_p=NULL;
72 17 : lowIm_p=NULL;
73 17 : lowImOrig_p=NULL;
74 17 : cwImage_p=NULL;
75 17 : cwHighIm_p=NULL;
76 17 : }
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 17 : Feather::~Feather(){
88 17 : highIm_p=NULL;
89 17 : lowIm_p=NULL;
90 17 : cwImage_p=NULL;
91 17 : }
92 :
93 17 : void Feather::setSDImage(const ImageInterface<Float>& SDImage){
94 34 : LogIO os(LogOrigin("Feather", "setSDImage()", WHERE));
95 17 : if(highIm_p.null())
96 0 : throw(AipsError("High res image has to be defined before SD image"));
97 17 : ImageInfo lowInfo=SDImage.imageInfo();
98 17 : lBeam_p=lowInfo.restoringBeam();
99 17 : if(lBeam_p.isNull())
100 0 : throw(AipsError("No Single dish restoring beam info in image"));
101 17 : CoordinateSystem csyslow=SDImage.coordinates();
102 17 : CountedPtr<ImageInterface<Float> > sdcopy;
103 17 : sdcopy=new TempImage<Float>(SDImage.shape(), csyslow);
104 17 : (*sdcopy).copyData(SDImage);
105 17 : ImageUtilities::copyMiscellaneous(*sdcopy, SDImage);
106 17 : if(SDImage.getDefaultMask() != "")
107 10 : Imager::copyMask(*sdcopy, SDImage, SDImage.getDefaultMask());
108 17 : std::unique_ptr<ImageInterface<Float> > copyPtr;
109 17 : std::unique_ptr<ImageInterface<Float> > copyPtr2;
110 :
111 :
112 17 : Vector<Stokes::StokesTypes> stokesvec;
113 17 : 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 17 : 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 17 : lowIm_p=new TempImage<Float>(highIm_p->shape(), csysHigh_p);
131 : // regrid the single dish image
132 : {
133 17 : Vector<Int> dirAxes=CoordinateUtil::findDirectionAxes(csysHigh_p);
134 17 : IPosition axes(3,dirAxes(0),dirAxes(1),2);
135 17 : Int spectralAxisIndex=CoordinateUtil::findSpectralAxis(csysHigh_p);
136 17 : axes(2)=spectralAxisIndex;
137 17 : if(sdcopy->getDefaultMask() != "")
138 10 : lowIm_p->makeMask(sdcopy->getDefaultMask(), true, true, true, true);
139 :
140 17 : ImageRegrid<Float> ir;
141 17 : ir.regrid(*lowIm_p, Interpolate2D::LINEAR, axes, *sdcopy);
142 17 : }
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 17 : if(lowImOrig_p.null()){
153 17 : lowImOrig_p=new TempImage<Float>(lowIm_p->shape(), lowIm_p->coordinates(),0);
154 17 : lowImOrig_p->copyData(*lowIm_p);
155 17 : lBeamOrig_p=lBeam_p;
156 : }
157 17 : cweightCalced_p=false;
158 17 : }
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 17 : void Feather::setINTImage(const ImageInterface<Float>& INTImage){
187 34 : LogIO os(LogOrigin("Feather", "setINTImage()", WHERE));
188 17 : ImageInfo highInfo=INTImage.imageInfo();
189 17 : hBeam_p=highInfo.restoringBeam();
190 17 : if(hBeam_p.isNull())
191 0 : throw(AipsError("No high-resolution restoring beam info in image"));
192 17 : csysHigh_p=INTImage.coordinates();
193 34 : CountedPtr<ImageInterface<Float> > intcopy=new TempImage<Float>(INTImage.shape(), csysHigh_p);
194 17 : (*intcopy).copyData(INTImage);
195 17 : ImageUtilities::copyMiscellaneous(*intcopy, INTImage);
196 17 : if(INTImage.getDefaultMask() != "")
197 0 : Imager::copyMask(*intcopy, INTImage, INTImage.getDefaultMask());
198 17 : std::unique_ptr<ImageInterface<Float> > copyPtr;
199 17 : std::unique_ptr<ImageInterface<Float> > copyPtr2;
200 :
201 :
202 17 : Vector<Stokes::StokesTypes> stokesvec;
203 17 : 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 17 : 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 17 : highIm_p=new TempImage<Float>(intcopy->shape(), csysHigh_p);
222 :
223 17 : highIm_p->copyData(*intcopy);
224 17 : ImageUtilities::copyMiscellaneous(*highIm_p, *intcopy);
225 17 : String maskname=intcopy->getDefaultMask();
226 17 : if(maskname != ""){
227 0 : Imager::copyMask(*highIm_p, *intcopy, maskname);
228 :
229 : }
230 17 : cweightCalced_p=false;
231 :
232 17 : }
233 :
234 3 : Bool Feather::setEffectiveDishDiam(const Float xdiam, const Float ydiam){
235 3 : dishDiam_p=ydiam >0 ? min(xdiam, ydiam) : xdiam;
236 : {//reset to the original SD image
237 3 : lowIm_p->copyData(*lowImOrig_p);
238 3 : lBeam_p=lBeamOrig_p;
239 : }
240 : //Change the beam of SD image
241 3 : Double freq = worldFreq(lowIm_p->coordinates(), 0);
242 : //cerr << "Freq " << freq << endl;
243 3 : Quantity halfpb=Quantity(1.22*C::c/freq/dishDiam_p, "rad");
244 6 : GaussianBeam newBeam(halfpb, halfpb, Quantity(0.0, "deg"));
245 3 : GaussianBeam toBeUsed;
246 : try {
247 3 : 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 2 : 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 2 : lBeam_p=newBeam;
261 : //reset cweight if it was calculated already
262 2 : cweightCalced_p=false;
263 2 : cweightApplied_p=false;
264 :
265 2 : return true;
266 5 : }
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 17 : void Feather::setSDScale(Float sdscale){
280 17 : sdScale_p=sdscale;
281 17 : }
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 16 : void Feather::applyFeather(){
375 16 : if(cweightApplied_p)
376 0 : return;
377 16 : calcCWeightImage();
378 16 : if(highIm_p.null())
379 0 : throw(AipsError("No high resolution image set"));
380 16 : cwHighIm_p=new TempImage<Complex>(highIm_p->shape(), highIm_p->coordinates() );
381 16 : StokesImageUtil::From(*cwHighIm_p, *highIm_p);
382 16 : LatticeFFT::cfft2d( *cwHighIm_p );
383 16 : Vector<Int> extraAxes(cwHighIm_p->shape().nelements()-2);
384 16 : if(extraAxes.nelements() > 0){
385 :
386 16 : if(extraAxes.nelements() ==2){
387 16 : Int n3=cwHighIm_p->shape()(2);
388 16 : Int n4=cwHighIm_p->shape()(3);
389 16 : IPosition blc(cwHighIm_p->shape());
390 16 : IPosition trc(cwHighIm_p->shape());
391 16 : blc(0)=0; blc(1)=0;
392 16 : trc(0)=cwHighIm_p->shape()(0)-1;
393 16 : trc(1)=cwHighIm_p->shape()(1)-1;
394 32 : for (Int j=0; j < n3; ++j){
395 34 : for (Int k=0; k < n4 ; ++k){
396 18 : blc(2)=j; trc(2)=j;
397 18 : blc(3)=k; trc(3)=k;
398 18 : Slicer sl(blc, trc, Slicer::endIsLast);
399 18 : SubImage<Complex> cimagehighSub(*cwHighIm_p, sl, true);
400 18 : cimagehighSub.copyData( (LatticeExpr<Complex>)((cimagehighSub * (*cwImage_p))));
401 18 : }
402 : }
403 16 : }
404 : }
405 : else{
406 0 : cwHighIm_p->copyData(
407 0 : (LatticeExpr<Complex>)(((*cwHighIm_p) * (*cwImage_p) )));
408 : }
409 16 : cweightApplied_p=true;
410 :
411 16 : }
412 :
413 16 : void Feather::calcCWeightImage(){
414 16 : if(cweightCalced_p)
415 0 : return;
416 16 : if(highIm_p.null())
417 0 : throw(AipsError("No Interferometer image was defined"));
418 16 : IPosition myshap(highIm_p->shape());
419 16 : Vector<Int> dirAxes;
420 16 : dirAxes=CoordinateUtil::findDirectionAxes(highIm_p->coordinates());
421 80 : for( uInt k=0; k< myshap.nelements(); ++k){
422 64 : if( (Int(k) != dirAxes[0]) && (Int(k) != dirAxes[1]))
423 32 : myshap(k)=1;
424 : }
425 :
426 16 : cwImage_p=new TempImage<Complex>(myshap, highIm_p->coordinates());
427 : {
428 16 : TempImage<Float> lowpsf(myshap, cwImage_p->coordinates());
429 16 : lowpsf.set(0.0);
430 32 : IPosition center(4, Int((myshap(0)/4)*2),
431 16 : Int((myshap(1)/4)*2),0,0);
432 16 : lowpsf.putAt(1.0, center);
433 16 : StokesImageUtil::Convolve(lowpsf, lBeam_p, false);
434 16 : StokesImageUtil::From(*cwImage_p, lowpsf);
435 16 : }
436 16 : LatticeFFT::cfft2d( *cwImage_p );
437 16 : LatticeExprNode node = max( *cwImage_p );
438 16 : Float fmax = abs(node.getComplex());
439 16 : cwImage_p->copyData( (LatticeExpr<Complex>)( 1.0f - (*cwImage_p)/fmax ) );
440 16 : cweightCalced_p=true;
441 16 : cweightApplied_p=false;
442 16 : }
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 0 : void Feather::saveFeatheredImage( const String& imagename ){
677 0 : PagedImage<Float> featherImage(highIm_p->shape(), highIm_p->coordinates(), imagename );
678 0 : saveFeatheredImage(imagename);
679 0 : }
680 :
681 :
682 16 : void Feather::saveFeatheredImage(ImageInterface<Float>& featherImage){
683 16 : applyFeather();
684 : Int stokesAxis, spectralAxis;
685 16 : spectralAxis=CoordinateUtil::findSpectralAxis(csysHigh_p);
686 16 : Vector<Stokes::StokesTypes> stokesvec;
687 16 : stokesAxis=CoordinateUtil::findStokesAxis(stokesvec, csysHigh_p);
688 16 : Vector<Int> dirAxes=CoordinateUtil::findDirectionAxes(csysHigh_p);
689 16 : Int n3=highIm_p->shape()(stokesAxis);
690 16 : Int n4=highIm_p->shape()(spectralAxis);
691 16 : IPosition blc(highIm_p->shape());
692 16 : IPosition trc(highIm_p->shape());
693 16 : blc(dirAxes(0))=0; blc(dirAxes(1))=0;
694 16 : trc(dirAxes(0))=highIm_p->shape()(dirAxes(0))-1;
695 16 : trc(dirAxes(1))=highIm_p->shape()(dirAxes(1))-1;
696 32 : TempImage<Complex>cimagelow(lowIm_p->shape(), lowIm_p->coordinates());
697 16 : StokesImageUtil::From(cimagelow, *lowIm_p);
698 16 : LatticeFFT::cfft2d( cimagelow);
699 16 : Float sdScaling = sdScale_p*hBeam_p.getArea("arcsec2")/lBeam_p.getArea("arcsec2");
700 32 : for (Int j=0; j < n3; ++j){
701 34 : for (Int k=0; k < n4 ; ++k){
702 18 : blc(stokesAxis)=j; trc(stokesAxis)=j;
703 18 : blc(spectralAxis)=k; trc(spectralAxis)=k;
704 18 : Slicer sl(blc, trc, Slicer::endIsLast);
705 18 : SubImage<Complex> cimagehighSub(*cwHighIm_p, sl, true);
706 18 : SubImage<Complex> cimagelowSub(cimagelow, sl, true);
707 18 : cimagelowSub.copyData( (LatticeExpr<Complex>)((cimagehighSub + cimagelowSub * sdScaling)));
708 18 : }
709 : }
710 : // FT back to image plane
711 16 : LatticeFFT::cfft2d( cimagelow, false);
712 :
713 : // write to output image
714 16 : StokesImageUtil::To(featherImage, cimagelow);
715 16 : ImageUtilities::copyMiscellaneous(featherImage, *highIm_p);
716 16 : String maskofHigh=highIm_p->getDefaultMask();
717 16 : String maskofLow=lowIm_p->getDefaultMask();
718 16 : if(maskofHigh != String("")){
719 : //ImageUtilities::copyMask(featherImage, *highIm_p, "mask0", maskofHigh, AxesSpecifier());
720 : //featherImage.setDefaultMask("mask0");
721 0 : Imager::copyMask(featherImage, *highIm_p, maskofHigh);
722 0 : if(maskofLow != String("")){
723 0 : featherImage.pixelMask().copyData(LatticeExpr<Bool> (featherImage.pixelMask() && (lowIm_p->pixelMask())));
724 :
725 : }
726 : }
727 16 : else if(maskofLow != String("")){
728 : //ImageUtilities::copyMask(featherImage, *lowIm_p, "mask0", maskofLow, AxesSpecifier());
729 9 : Imager::copyMask(featherImage, *lowIm_p, maskofLow);
730 : }
731 16 : }
732 :
733 0 : void Feather::applyDishDiam(ImageInterface<Complex>& image, GaussianBeam& beam, Float effDiam, ImageInterface<Float>& fftim, Vector<Quantity>& extraconv){
734 : /*
735 : MathFunc<Float> dd(SPHEROIDAL);
736 : Vector<Float> valsph(31);
737 : for(Int k=0; k <31; ++k){
738 : valsph(k)=dd.value((Float)(k)/10.0);
739 : }
740 : Quantity fulrad((52.3/2.0*3.0/1.1117),"arcsec");
741 : Quantity lefreq(224.0/effDiam*9.0, "GHz");
742 :
743 : PBMath1DNumeric elpb(valsph, fulrad, lefreq);*/
744 : ////////////////////////
745 : /*{
746 : PagedImage<Float> thisScreen(image.shape(), image.coordinates(), "Before_apply.image");
747 : LatticeExpr<Float> le(abs(image));
748 : thisScreen.copyData(le);
749 : }*/
750 : //////////////////////
751 :
752 0 : Double freq = worldFreq(image.coordinates(), 0);
753 : //cerr << "Freq " << freq << endl;
754 0 : Quantity halfpb=Quantity(1.22*C::c/freq/effDiam, "rad");
755 : //PBMath1DAiry elpb( Quantity(effDiam,"m"), Quantity(0.01,"m"),
756 : // Quantity(0.8564,"deg"), Quantity(1.0,"GHz"));
757 0 : GaussianBeam newBeam(halfpb, halfpb, Quantity(0.0, "deg"));
758 : try {
759 :
760 0 : GaussianBeam toBeUsed;
761 : //cerr << "beam " << beam.toVector() << endl;
762 : //cerr << "newBeam " << newBeam.toVector() << endl;
763 0 : GaussianDeconvolver::deconvolve(toBeUsed, newBeam, beam);
764 0 : extraconv.resize(3);
765 : // use the Major difference
766 0 : extraconv(0) = toBeUsed.getMajor();
767 0 : extraconv(1) = toBeUsed.getMinor();
768 0 : extraconv(2) = toBeUsed.getPA();
769 :
770 0 : }
771 0 : catch (const AipsError& x) {
772 0 : throw(AipsError("Beam due to new effective diameter may be smaller than the beam of original dish image"));
773 0 : }
774 : //////////////////////
775 : //1 GHz equivalent
776 : // halfpb=Quantity(1.22*C::c/1.0e9/effDiam, "rad");
777 : //cerr << "halfpb " << halfpb << endl;
778 : //PBMath1DGauss elpb(halfpb, Quantity(0.8564,"deg"), Quantity(1.0,"GHz"), true);
779 :
780 0 : fftim.set(0.0);
781 :
782 0 : IPosition center(4, Int((fftim.shape()(0)/4)*2),
783 0 : Int((fftim.shape()(1)/4)*2),0,0);
784 0 : fftim.putAt(1.0, center);
785 0 : StokesImageUtil::Convolve(fftim, newBeam, false);
786 0 : StokesImageUtil::From(image, fftim);
787 : /*
788 : TempImage<Complex> elbeamo(image.shape(), image.coordinates());
789 : elbeamo.set(1.0);
790 : MDirection wcenter;
791 : {
792 : Int directionIndex=
793 : image.coordinates().findCoordinate(Coordinate::DIRECTION);
794 : DirectionCoordinate
795 : directionCoord=image.coordinates().directionCoordinate(directionIndex);
796 : Vector<Double> pcenter(2);
797 : pcenter(0) = image.shape()(directionIndex)/2;
798 : pcenter(1) = image.shape()(directionIndex+1)/2;
799 : directionCoord.toWorld( wcenter, pcenter );
800 : }
801 : elpb.applyPB(elbeamo, elbeamo, wcenter, Quantity(0.0, "deg"),
802 : BeamSquint::NONE);
803 :
804 :
805 : StokesImageUtil::To(fftim, elbeamo);
806 : */
807 : ///////////////////
808 : //StokesImageUtil::FitGaussianPSF(fftim,
809 : // beam);
810 : //cerr << "To be convolved beam 2 " << beam.toVector() << endl;
811 : ////////////////
812 :
813 0 : LatticeFFT::cfft2d(image);
814 : //image.copyData((LatticeExpr<Complex>)(elbeamo) );
815 : //elbeamo.copyData(image);
816 : // LatticeFFT::cfft2d(elbeamo, false);
817 : //StokesImageUtil::To(fftim, elbeamo);
818 : //StokesImageUtil::FitGaussianPSF(fftim,
819 : // beam);
820 : /////////
821 : /*{
822 : PagedImage<Float> thisScreen(image.shape(), image.coordinates(), "After_apply.image");
823 : //LatticeExpr<Float> le(abs(image));
824 : thisScreen.copyData(fftim);
825 : }*/
826 : ////////
827 0 : }
828 :
829 :
830 :
831 17 : void Feather::feather(ImageInterface<Float>& featherImage, 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){
832 :
833 17 : Feather plume;
834 17 : plume.setINTImage(high);
835 17 : ImageInfo lowInfo=low0.imageInfo();
836 17 : GaussianBeam lBeam=lowInfo.restoringBeam();
837 17 : if(lBeam.isNull()) {
838 0 : getLowBeam(low0, lowPSF, useDefaultPB, vpTableStr, lBeam);
839 0 : if(lBeam.isNull())
840 0 : throw(AipsError("Do not have low resolution beam info "));
841 0 : TempImage<Float> newLow(low0.shape(), low0.coordinates());
842 0 : newLow.copyData(low0);
843 0 : lowInfo.removeRestoringBeam();
844 0 : lowInfo.setRestoringBeam(lBeam);
845 0 : newLow.setImageInfo(lowInfo);
846 0 : plume.setSDImage(newLow);
847 0 : }
848 : else{
849 17 : plume.setSDImage(low0);
850 : }
851 17 : plume.setSDScale(sdScale);
852 17 : if(effDiam >0.0){
853 2 : plume.setEffectiveDishDiam(effDiam, effDiam);
854 : }
855 15 : else if(doHiPassFilterOnSD){
856 : Float xdiam, ydiam;
857 1 : plume.getEffectiveDishDiam(xdiam, ydiam);
858 1 : plume.setEffectiveDishDiam(xdiam, ydiam);
859 : }
860 16 : plume.saveFeatheredImage(featherImage);
861 :
862 19 : }
863 :
864 :
865 0 : void Feather::getLowBeam(const ImageInterface<Float>& low0, const String& lowPSF, const Bool useDefaultPB, const String& vpTableStr, GaussianBeam& lBeam){
866 0 : LogIO os(LogOrigin("Feather", "getLowBeam()", WHERE));
867 0 : PBMath * myPBp = 0;
868 0 : if((lowPSF=="") && lBeam.isNull()) {
869 : // create the low res's PBMath object, needed to apply PB
870 : // to make high res Fourier weight image
871 0 : if (useDefaultPB) {
872 : // look up the telescope in ObsInfo
873 0 : ObsInfo oi = low0.coordinates().obsInfo();
874 0 : String myTelescope = oi.telescope();
875 0 : if (myTelescope == "") {
876 : os << LogIO::SEVERE << "No telescope imbedded in low res image"
877 0 : << LogIO::POST;
878 : os << LogIO::SEVERE << "Create a PB description with the vpmanager"
879 0 : << LogIO::EXCEPTION;
880 : }
881 0 : Quantity qFreq;
882 : {
883 0 : Double freq = worldFreq(low0.coordinates(), 0);
884 0 : qFreq = Quantity( freq, "Hz" );
885 : }
886 0 : String band;
887 : PBMath::CommonPB whichPB;
888 0 : String pbName;
889 : // get freq from coordinates
890 0 : PBMath::whichCommonPBtoUse (myTelescope, qFreq, band, whichPB,
891 : pbName);
892 0 : if (whichPB == PBMath::UNKNOWN) {
893 : os << LogIO::SEVERE << "Unknown telescope for PB type: "
894 0 : << myTelescope << LogIO::EXCEPTION;
895 : }
896 0 : myPBp = new PBMath(whichPB);
897 0 : } else {
898 : // get the PB from the vpTable
899 0 : Table vpTable( vpTableStr );
900 0 : ScalarColumn<TableRecord> recCol(vpTable, (String)"pbdescription");
901 0 : myPBp = new PBMath(recCol(0));
902 0 : }
903 0 : AlwaysAssert((myPBp != 0), AipsError);
904 : }
905 :
906 :
907 : // get image center direction
908 0 : MDirection wcenter;
909 : {
910 : Int directionIndex=
911 0 : low0.coordinates().findCoordinate(Coordinate::DIRECTION);
912 0 : AlwaysAssert(directionIndex>=0, AipsError);
913 : DirectionCoordinate
914 0 : directionCoord=low0.coordinates().directionCoordinate(directionIndex);
915 0 : Vector<Double> pcenter(2);
916 0 : pcenter(0) = low0.shape()(0)/2;
917 0 : pcenter(1) = low0.shape()(1)/2;
918 0 : directionCoord.toWorld( wcenter, pcenter );
919 0 : }
920 :
921 : // make the weight image
922 0 : IPosition myshap(low0.shape());
923 0 : for( uInt k=2; k< myshap.nelements(); ++k){
924 0 : myshap(k)=1;
925 : }
926 :
927 0 : TempImage<Float> lowpsf0;
928 0 : TempImage<Complex> cweight(myshap, low0.coordinates());
929 0 : if(lowPSF=="") {
930 : os << LogIO::NORMAL // Loglevel INFO
931 0 : << "Using primary beam to determine weighting.\n" << LogIO::POST;
932 0 : cweight.set(1.0);
933 0 : if (myPBp != 0) {
934 0 : myPBp->applyPB(cweight, cweight, wcenter, Quantity(0.0, "deg"),
935 : BeamSquint::NONE);
936 :
937 0 : lowpsf0=TempImage<Float>(cweight.shape(), cweight.coordinates());
938 :
939 : os << LogIO::NORMAL // Loglevel INFO
940 : << "Determining scaling from SD Primary Beam.\n"
941 0 : << LogIO::POST;
942 0 : StokesImageUtil::To(lowpsf0, cweight);
943 0 : StokesImageUtil::FitGaussianPSF(lowpsf0,
944 : lBeam);
945 : }
946 0 : delete myPBp;
947 : }
948 : else {
949 : os << LogIO::NORMAL // Loglevel INFO
950 : << "Using specified low resolution PSF to determine weighting.\n"
951 0 : << LogIO::POST;
952 : // regrid the single dish psf
953 0 : PagedImage<Float> lowpsfDisk(lowPSF);
954 0 : IPosition lshape(lowpsfDisk.shape());
955 0 : lshape.resize(4);
956 0 : lshape(2)=1; lshape(3)=1;
957 0 : TempImage<Float>lowpsf(lshape,lowpsfDisk.coordinates());
958 0 : IPosition blc(lowpsfDisk.shape());
959 0 : IPosition trc(lowpsfDisk.shape());
960 0 : blc(0)=0; blc(1)=0;
961 0 : trc(0)=lowpsfDisk.shape()(0)-1;
962 0 : trc(1)=lowpsfDisk.shape()(1)-1;
963 0 : for( uInt k=2; k < lowpsfDisk.shape().nelements(); ++k){
964 0 : blc(k)=0; trc(k)=0;
965 : }// taking first plane
966 0 : Slicer sl(blc, trc, Slicer::endIsLast);
967 0 : lowpsf.copyData(SubImage<Float>(lowpsfDisk, sl, false));
968 : os << LogIO::NORMAL // Loglevel INFO
969 0 : << "Determining scaling from low resolution PSF.\n" << LogIO::POST;
970 0 : StokesImageUtil::FitGaussianPSF(lowpsf, lBeam);
971 0 : }
972 :
973 :
974 :
975 0 : }
976 :
977 4 : Double Feather::worldFreq(const CoordinateSystem& cs, Int spectralpix){
978 : ///freq part
979 4 : Int spectralIndex=cs.findCoordinate(Coordinate::SPECTRAL);
980 : SpectralCoordinate
981 : spectralCoord=
982 4 : cs.spectralCoordinate(spectralIndex);
983 4 : Vector<String> units(1); units = "Hz";
984 4 : spectralCoord.setWorldAxisUnits(units);
985 4 : Vector<Double> spectralWorld(1);
986 4 : Vector<Double> spectralPixel(1);
987 4 : spectralPixel(0) = spectralpix;
988 4 : spectralCoord.toWorld(spectralWorld, spectralPixel);
989 4 : Double freq = spectralWorld(0);
990 4 : return freq;
991 4 : }
992 :
993 :
994 :
995 :
996 :
997 : }//# NAMESPACE CASA - END
|
