Line data Source code
1 : // -*- C++ -*-
2 : //# WOnlyConvFunc.cc: Implementation of the WOnlyConvFunc class
3 : //# Copyright (C) 1997,1998,1999,2000,2001,2002,2003
4 : //# Associated Universities, Inc. Washington DC, USA.
5 : //#
6 : //# This library is free software; you can redistribute it and/or modify it
7 : //# under the terms of the GNU Library General Public License as published by
8 : //# the Free Software Foundation; either version 2 of the License, or (at your
9 : //# option) any later version.
10 : //#
11 : //# This library is distributed in the hope that it will be useful, but WITHOUT
12 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
13 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
14 : //# License for more details.
15 : //#
16 : //# You should have received a copy of the GNU Library General Public License
17 : //# along with this library; if not, write to the Free Software Foundation,
18 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
19 : //#
20 : //# Correspondence concerning AIPS++ should be addressed as follows:
21 : //# Internet email: casa-feedback@nrao.edu.
22 : //# Postal address: AIPS++ Project Office
23 : //# National Radio Astronomy Observatory
24 : //# 520 Edgemont Road
25 : //# Charlottesville, VA 22903-2475 USA
26 : //#
27 : //# $Id$
28 : //
29 : #include <synthesis/TransformMachines/WOnlyConvFunc.h>
30 : #include <synthesis/TransformMachines/SynthesisError.h>
31 : #include <synthesis/TransformMachines/WTerm.h>
32 : #include <synthesis/TransformMachines/PSTerm.h>
33 : #include <casacore/images/Images/ImageInterface.h>
34 : #include <synthesis/TransformMachines/Utils.h>
35 : #include <synthesis/TransformMachines/CFStore.h>
36 : #include <synthesis/TransformMachines/ConvolutionFunction.h>
37 : #include <casacore/coordinates/Coordinates/DirectionCoordinate.h>
38 : #include <casacore/coordinates/Coordinates/SpectralCoordinate.h>
39 : #include <casacore/coordinates/Coordinates/StokesCoordinate.h>
40 : #include <casacore/lattices/LatticeMath/LatticeFFT.h>
41 : #include <casacore/casa/ostream.h>
42 : using namespace casacore;
43 : namespace casa{
44 0 : void WOnlyConvFunc::makeConvFunction(const ImageInterface<Complex>& image,
45 : const VisBuffer& vb,
46 : const Int wConvSize,
47 : const Float pa,
48 : CFStore& cfs,
49 : CFStore& cfwts)
50 : {
51 0 : LogIO log_l(LogOrigin("WOnlyConvFunc", "makeConvFunction"));
52 0 : Int convSize=0, convSampling=4;
53 0 : Double wScale=1;
54 0 : Array<Complex> convFunc_l;
55 0 : Double cfRefFreq=-1.0;
56 :
57 0 : Int nx=image.shape()(0);
58 :
59 0 : if(wConvSize>0)
60 : {
61 0 : log_l << "Using " << wConvSize << " planes for W-projection" << LogIO::POST;
62 : Double maxUVW;
63 0 : maxUVW=0.25/abs(image.coordinates().increment()(0));
64 : log_l << "Estimating maximum possible W = " << maxUVW
65 0 : << " (wavelengths)" << LogIO::POST;
66 :
67 0 : Double invLambdaC=vb.frequency()(0)/C::c;
68 0 : Double invMinL = vb.frequency()((vb.frequency().nelements())-1)/C::c;
69 : log_l << "wavelength range = " << 1.0/invLambdaC << " (m) to "
70 0 : << 1.0/invMinL << " (m)" << LogIO::POST;
71 0 : if (wConvSize > 1)
72 : {
73 0 : wScale=Float((wConvSize-1)*(wConvSize-1))/maxUVW;
74 : log_l << "Scaling in W (at maximum W) = " << 1.0/wScale
75 0 : << " wavelengths per pixel" << LogIO::POST;
76 : }
77 : }
78 : // Get the coordinate system
79 : //
80 0 : CoordinateSystem coords(image.coordinates());
81 : //
82 : // Make a two dimensional image to calculate auto-correlation of
83 : // the ideal illumination pattern. We want this on a fine grid in
84 : // the UV plane
85 : //
86 0 : Int directionIndex=coords.findCoordinate(Coordinate::DIRECTION);
87 0 : AlwaysAssert(directionIndex>=0, AipsError);
88 0 : DirectionCoordinate dc=coords.directionCoordinate(directionIndex);
89 0 : Vector<Double> sampling;
90 0 : sampling = dc.increment();
91 :
92 0 : sampling*=Double(convSampling);
93 0 : sampling*=Double(nx)/Double(convSize);
94 :
95 0 : dc.setIncrement(sampling);
96 :
97 0 : Vector<Double> unitVec(2);
98 0 : unitVec=convSize/2;
99 0 : dc.setReferencePixel(unitVec);
100 :
101 : // Set the reference value to that of the image
102 0 : coords.replaceCoordinate(dc, directionIndex);
103 0 : IPosition pbShape(4, convSize, convSize, 1, 1);
104 0 : TempImage<Complex> twoDPB(pbShape, coords);
105 :
106 0 : Int inner=convSize/convSampling, polInUse = 1;
107 0 : cfs.data = new Array<Complex>(IPosition(4,convSize,convSize, wConvSize,polInUse));
108 0 : convFunc_l.reference(*cfs.data); convFunc_l=0;
109 0 : cfs.resize(wConvSize);
110 :
111 0 : IPosition start(4, 0, 0, 0, 0);
112 0 : IPosition pbSlice(4, convSize, convSize, 1, 1);
113 :
114 0 : Matrix<Complex> screen(convSize, convSize);
115 :
116 0 : for (Int iw=0;iw<wConvSize;iw++)
117 : {
118 0 : wTerm_p->applySky(screen, iw, sampling, wScale, convSize);
119 :
120 0 : psTerm_p->applySky(screen, /*iw,*/ sampling, /*wScale,*/ inner);
121 :
122 : //
123 : // Fill the complex image with the w-term...
124 : //
125 0 : IPosition PolnPlane(4,0,0,0,0);
126 0 : IPosition ndx(4,0,0,0,0);
127 :
128 0 : for(Int i=0;i<polInUse;i++)
129 : {
130 0 : PolnPlane(2)=i;
131 0 : twoDPB.putSlice(screen, PolnPlane);
132 : }
133 :
134 0 : Complex cpeak=max(twoDPB.get());
135 0 : twoDPB.put(twoDPB.get()/cpeak);
136 :
137 0 : CoordinateSystem cs=twoDPB.coordinates();
138 0 : Int index= twoDPB.coordinates().findCoordinate(Coordinate::SPECTRAL);
139 0 : SpectralCoordinate SpCS = twoDPB.coordinates().spectralCoordinate(index);
140 :
141 0 : cfRefFreq=SpCS.referenceValue()(0);
142 0 : Vector<Double> refValue; refValue.resize(1); refValue(0)=cfRefFreq;
143 0 : SpCS.setReferenceValue(refValue);
144 0 : cs.replaceCoordinate(SpCS,index);
145 : //
146 : // Now FFT and get the result back
147 : //
148 0 : LatticeFFT::cfft2d(twoDPB);
149 : //
150 : // Fill the convolution function planes with the result.
151 : //
152 : {
153 0 : IPosition start(4, 0, 0, 0, 0),
154 0 : pbSlice(4, twoDPB.shape()[0]-1, twoDPB.shape()[1]-1, polInUse, 1);
155 0 : IPosition sliceStart(4,0,0,iw,0),
156 0 : sliceLength(4,convFunc_l.shape()[0]-1,convFunc_l.shape()[1]-1,1,polInUse);
157 :
158 0 : convFunc_l(Slicer(sliceStart,sliceLength)).nonDegenerate()
159 0 : =(twoDPB.getSlice(start, pbSlice, true));
160 0 : }
161 0 : }
162 : //
163 : // Set various CF parameters
164 : //
165 0 : cfs.sampling(0)=convSampling;
166 0 : setSupport(convFunc_l, cfs);
167 :
168 0 : Int index=coords.findCoordinate(Coordinate::SPECTRAL);
169 0 : SpectralCoordinate spCS = coords.spectralCoordinate(index);
170 0 : Vector<Double> refValue; refValue.resize(1);refValue(0)=cfRefFreq;
171 0 : spCS.setReferenceValue(refValue);
172 0 : coords.replaceCoordinate(spCS,index);
173 :
174 0 : cfs.coordSys=coords;
175 :
176 : (void)pa; (void)cfwts; // To suppress compiler warnings
177 0 : };
178 :
179 0 : void WOnlyConvFunc::setSupport(Array<Complex>& convFunc, CFStore& cfs)
180 : {
181 0 : Int maxConvSupport=-1, polInUse = 1;
182 0 : Int ConvFuncOrigin=convFunc.shape()[0]/2; // Conv. Func. is half that size of convSize
183 0 : IPosition ndx(4,ConvFuncOrigin,0,0,0), cfShape(convFunc.shape());
184 :
185 : //UNUSED: Int maxConvWtSupport=0, supportBuffer;
186 0 : for (Int iw=0;iw<cfShape[2];iw++)
187 : {
188 0 : Bool found=false;
189 : Float threshold;
190 : Int R;
191 0 : ndx(2) = iw;
192 :
193 0 : ndx(0)=ndx(1)=ConvFuncOrigin;
194 0 : ndx(2) = iw;
195 0 : threshold = abs(convFunc(ndx))*wTerm_p->getSupportThreshold();
196 : //
197 : // Find the support size of the conv. function in pixels
198 : //
199 : //UNUSED: Int wtR;
200 0 : found = findSupport(convFunc,threshold,ConvFuncOrigin,R);
201 : //
202 : // Set the support size for each W-plane and for all
203 : // Pol-planes. Assuming that support size for all Pol-planes
204 : // is same.
205 : //
206 0 : if(found)
207 : {
208 : // Int maxR=R;//max(ndx(0),ndx(1));
209 0 : for(Int ipol=0;ipol<polInUse;ipol++)
210 : {
211 0 : cfs.xSupport(iw)=cfs.ySupport(iw)=Int(R/cfs.sampling(0));
212 0 : cfs.xSupport(iw)=cfs.ySupport(iw)=Int(0.5+Float(R)/cfs.sampling(0))+1;
213 :
214 0 : if (cfs.maxXSupport == -1)
215 0 : if (cfs.xSupport(iw) > maxConvSupport)
216 0 : maxConvSupport = cfs.xSupport(iw);
217 : }
218 : }
219 : }
220 0 : }
221 :
222 : //
223 : //-----------------------------------------------------------------------
224 : //
225 0 : Bool WOnlyConvFunc::findSupport(Array<Complex>& func, Float& threshold,
226 : Int& origin, Int& R)
227 : {
228 0 : LogIO log_l(LogOrigin("WOnlyConvFunc", "findSupport"));
229 : Double NSteps;
230 0 : Int PixInc=1;
231 0 : Vector<Complex> vals;
232 0 : IPosition ndx(4,origin,0,0,0);
233 0 : Bool found=false;
234 0 : IPosition cfShape=func.shape();
235 0 : Int convSize = cfShape(0);
236 0 : for(R=convSize/4;R>1;R--)
237 : {
238 0 : NSteps = 90*R/PixInc; //Check every PixInc pixel along a
239 : //circle of radious R
240 0 : vals.resize((Int)(NSteps+0.5));
241 0 : vals=0;
242 0 : for(Int th=0;th<NSteps;th++)
243 : {
244 0 : ndx(0)=(int)(origin + R*sin(2.0*M_PI*th*PixInc/R));
245 0 : ndx(1)=(int)(origin + R*cos(2.0*M_PI*th*PixInc/R));
246 :
247 0 : if ((ndx(0) < cfShape(0)) && (ndx(1) < cfShape(1)))
248 0 : vals(th)=func(ndx);
249 : }
250 0 : if (max(abs(vals)) > threshold)
251 0 : {found=true;break;}
252 : }
253 0 : return found;
254 0 : };
255 :
256 : };
|
