Line data Source code
1 : //# WPConvFunc.cc: implementation of WPConvFunc
2 : //# Copyright (C) 2007
3 : //# Associated Universities, Inc. Washington DC, USA.
4 : //#
5 : //# This library is free software; you can redistribute it and/or modify it
6 : //# under the terms of the GNU Library General Public License as published by
7 : //# the Free Software Foundation; either version 2 of the License, or (at your
8 : //# option) any later version.
9 : //#
10 : //# This library 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 Library General Public
13 : //# License for more details.
14 : //#
15 : //# You should have received a copy of the GNU Library General Public License
16 : //# along with this library; if not, write to the Free Software Foundation,
17 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
18 : //#
19 : //# Correspondence concerning AIPS++ should be adressed as follows:
20 : //# Internet email: aips2-request@nrao.edu.
21 : //# Postal address: AIPS++ Project Office
22 : //# National Radio Astronomy Observatory
23 : //# 520 Edgemont Road
24 : //# Charlottesville, VA 22903-2475 USA
25 : //#
26 : //#
27 : //# $Id$
28 : #include <sstream>
29 : #include <iostream>
30 : #include <iomanip>
31 : #include <casacore/casa/Arrays/ArrayMath.h>
32 : #include <casacore/casa/Arrays/Array.h>
33 : #include <casacore/casa/Arrays/MaskedArray.h>
34 : #include <casacore/casa/Arrays/Vector.h>
35 : #include <casacore/casa/Arrays/Slice.h>
36 : #include <casacore/casa/Arrays/Matrix.h>
37 : #include <casacore/casa/Arrays/Cube.h>
38 : #include <casacore/casa/OS/HostInfo.h>
39 : #include <casacore/casa/OS/Timer.h>
40 : #include <casacore/casa/Utilities/Assert.h>
41 : #include <casacore/casa/Utilities/CompositeNumber.h>
42 : #include <casacore/coordinates/Coordinates/CoordinateSystem.h>
43 : #include <casacore/coordinates/Coordinates/DirectionCoordinate.h>
44 :
45 : #include <casacore/images/Images/ImageInterface.h>
46 : #include <casacore/images/Images/PagedImage.h>
47 : #include <casacore/images/Images/TempImage.h>
48 : #include <casacore/casa/Logging/LogIO.h>
49 : #include <casacore/casa/Logging/LogSink.h>
50 : #include <casacore/casa/Logging/LogMessage.h>
51 :
52 : #include <casacore/lattices/Lattices/ArrayLattice.h>
53 : #include <casacore/lattices/Lattices/SubLattice.h>
54 : #include <casacore/lattices/LRegions/LCBox.h>
55 : #include <casacore/lattices/LEL/LatticeExpr.h>
56 : #include <casacore/lattices/Lattices/LatticeCache.h>
57 : #include <casacore/lattices/LatticeMath/LatticeFFT.h>
58 : #include <casacore/scimath/Mathematics/ConvolveGridder.h>
59 : #include <casacore/scimath/Mathematics/FFTPack.h>
60 : #include <msvis/MSVis/VisBuffer.h>
61 : #include <msvis/MSVis/VisibilityIterator.h>
62 : #include <synthesis/TransformMachines/SimplePBConvFunc.h> //por SINCOS
63 :
64 : #include <synthesis/TransformMachines/WPConvFunc.h>
65 : #ifdef _OPENMP
66 : #include <omp.h>
67 : #endif
68 :
69 :
70 : using namespace casacore;
71 : namespace casa { //# NAMESPACE CASA - BEGIN
72 :
73 : typedef unsigned long long ooLong;
74 :
75 0 : WPConvFunc::WPConvFunc(const Double minW, const Double maxW, const Double rmsW):
76 0 : actualConvIndex_p(-1), convSize_p(0), convSupport_p(0), minW_p(minW), maxW_p(maxW), rmsW_p(rmsW) {
77 : //
78 0 : }
79 :
80 0 : WPConvFunc::~WPConvFunc(){
81 : //usage of CountedPtr keeps this simple
82 :
83 0 : }
84 :
85 0 : WPConvFunc::WPConvFunc(const RecordInterface& rec):
86 0 : actualConvIndex_p(-1), convSize_p(0), convSupport_p(0){
87 0 : String error;
88 0 : if (!fromRecord(error, rec)) {
89 0 : throw (AipsError("Failed to create WPConvFunc: " + error));
90 : }
91 :
92 0 : }
93 0 : WPConvFunc::WPConvFunc(const WPConvFunc& other):
94 0 : actualConvIndex_p(-1), convSize_p(0), convSupport_p(0){
95 :
96 0 : operator=(other);
97 0 : }
98 :
99 0 : WPConvFunc& WPConvFunc::operator=(const WPConvFunc& other){
100 0 : if(this != &other){
101 0 : uInt numConv=other.convFunctions_p.nelements();
102 0 : convFunctions_p.resize(numConv, true, false);
103 0 : convSupportBlock_p.resize(numConv, true, false);
104 0 : for (uInt k=0; k < numConv; ++k){
105 0 : convFunctions_p[k]=new Cube<Complex>();
106 0 : *(convFunctions_p[k])=*(other.convFunctions_p[k]);
107 0 : convSupportBlock_p[k]=new Vector<Int> ();
108 0 : *(convSupportBlock_p[k]) = *(other.convSupportBlock_p[k]);
109 : }
110 :
111 0 : convFunctionMap_p=other.convFunctionMap_p;
112 0 : convSizes_p.resize();
113 0 : convSizes_p=other.convSizes_p;
114 0 : actualConvIndex_p=other.actualConvIndex_p;
115 0 : convSize_p=other.convSize_p;
116 0 : convSupport_p.resize();
117 0 : convSupport_p=other.convSupport_p;
118 0 : convFunc_p.resize();
119 0 : convFunc_p=other.convFunc_p;
120 0 : wScaler_p=other.wScaler_p;
121 0 : convSampling_p=other.convSampling_p;
122 0 : nx_p=other.nx_p;
123 0 : ny_p=other.ny_p;
124 0 : minW_p=other.minW_p;
125 0 : maxW_p=other.maxW_p;
126 0 : rmsW_p=other.rmsW_p;
127 :
128 :
129 :
130 : }
131 0 : return *this;
132 : }
133 :
134 0 : void WPConvFunc::findConvFunction(const ImageInterface<Complex>& image,
135 : const VisBuffer& vb,
136 : const Int& wConvSizeUser,
137 : const Vector<Double>& uvScale,
138 : const Vector<Double>& uvOffset,
139 : const Float& padding,
140 : Int& convSampling,
141 : Cube<Complex>& convFunc,
142 : Int& convSize,
143 : Vector<Int>& convSupport, Double& wScale){
144 0 : if(checkCenterPix(image)){
145 0 : convFunc.resize();
146 0 : convFunc.reference(convFunc_p);
147 0 : convSize=convSize_p;
148 0 : convSampling=convSampling_p;
149 0 : convSupport.resize();
150 0 : convSupport=convSupport_p;
151 0 : wScale=Float((convFunc.shape()(2)-1)*(convFunc.shape()(2)-1))/wScaler_p;
152 0 : return;
153 : }
154 0 : LogIO os;
155 0 : os << LogOrigin("WPConvFunc", "findConvFunction") << LogIO::NORMAL;
156 :
157 :
158 : // Get the coordinate system
159 0 : CoordinateSystem coords(image.coordinates());
160 0 : Int directionIndex=coords.findCoordinate(Coordinate::DIRECTION);
161 0 : nx_p=Int(image.shape()(directionIndex));
162 0 : ny_p=Int(image.shape()(directionIndex+1));
163 :
164 0 : Int wConvSize=wConvSizeUser;
165 : ////Automatic mode
166 : Double maxUVW;
167 0 : if(wConvSize < 1){
168 0 : maxUVW=rmsW_p < 0.5*(minW_p+maxW_p) ? 1.05*maxW_p: (rmsW_p /(0.5*((minW_p)+maxW_p))*1.05*maxW_p) ;
169 0 : wConvSize=Int(maxUVW*fabs(sin(fabs(image.coordinates().increment()(0))*max(nx_p, ny_p)/2.0)));
170 0 : convSupport.resize(wConvSize);
171 : }
172 : else{
173 0 : if(maxW_p> 0.0)
174 0 : maxUVW= 1.05*maxW_p;
175 : else
176 0 : maxUVW=0.25/abs(image.coordinates().increment()(0));
177 :
178 : }
179 0 : if(wConvSize>1) {
180 0 : os << "W projection using " << wConvSize << " planes" << LogIO::POST;
181 :
182 : os << "Using maximum possible W = " << maxUVW
183 0 : << " (wavelengths)" << LogIO::POST;
184 :
185 0 : Double invLambdaC=vb.frequency()(0)/C::c;
186 : os << "Typical wavelength = " << 1.0/invLambdaC
187 0 : << " (m)" << LogIO::POST;
188 :
189 :
190 0 : wScale=Float((wConvSize-1)*(wConvSize-1))/maxUVW;
191 : //wScale=Float(wConvSize-1)/maxUVW;
192 0 : wScaler_p=maxUVW;;
193 0 : os << "Scaling in W (at maximum W) = " << 1.0/wScale
194 0 : << " wavelengths per pixel" << LogIO::POST;
195 : }
196 :
197 0 : Int planesPerChunk=100;
198 :
199 :
200 0 : if(wConvSize>1) {
201 0 : Int maxMemoryMB=min(uInt64(HostInfo::memoryTotal(true)), uInt64(HostInfo::memoryFree()))/1024;
202 : ////Common you have to have 4 GB
203 0 : if(maxMemoryMB < 4000)
204 0 : maxMemoryMB=4000;
205 0 : convSize=max(Int(nx_p*padding),Int(ny_p*padding));
206 : ///Do the same thing as in WProject::init
207 0 : CompositeNumber cn(convSize);
208 0 : convSize = cn.nextLargerEven(convSize);
209 : //nominal 100 wprojplanes above that you may (or not) go swapping
210 : //cerr << "maxmem " << maxMemoryMB << " npixels " << sqrt(Double(maxMemoryMB)/8.0*1024.0*1024.0) << endl;
211 0 : planesPerChunk=Int((Double(maxMemoryMB)/8.0*1024.0*1024.0)/Double(convSize)/Double(convSize));
212 : //cerr << "planesPerChunk " << planesPerChunk << endl;
213 0 : planesPerChunk=min(planesPerChunk, 100);
214 : // CompositeNumber cn(Int(maxConvSizeConsidered/2.0)*2);
215 0 : convSampling_p=4;
216 :
217 : //convSize=min(convSize,(Int)cn.nearestEven(Int(maxConvSizeConsidered/2.0)*2));
218 :
219 : }
220 : else{
221 0 : convSampling_p=1;
222 0 : convSize=max(Int(nx_p*padding),Int(ny_p*padding));
223 : }
224 0 : convSampling=convSampling_p;
225 0 : Int maxConvSize=convSize;
226 0 : convSupport.resize(wConvSize);
227 0 : convSupport.set(-1);
228 : ////set sampling
229 0 : AlwaysAssert(directionIndex>=0, AipsError);
230 0 : DirectionCoordinate dc=coords.directionCoordinate(directionIndex);
231 0 : Vector<Double> sampling;
232 0 : sampling = dc.increment();
233 0 : sampling*=Double(convSampling_p);
234 : //sampling*=Double(max(nx,ny))/Double(convSize);
235 0 : sampling[0]*=Double(padding*Float(nx_p))/Double(convSize);
236 0 : sampling[1]*=Double(padding*Float(ny_p))/Double(convSize);
237 : /////
238 0 : Int inner=convSize/convSampling_p;
239 : ConvolveGridder<Double, Complex>
240 0 : ggridder(IPosition(2, inner, inner), uvScale, uvOffset, "SF");
241 :
242 0 : Int nchunk=wConvSize/planesPerChunk;
243 0 : if((wConvSize%planesPerChunk) !=0)
244 0 : nchunk +=1;
245 0 : Vector<Int> chunksize(nchunk, planesPerChunk);
246 0 : if((wConvSize%planesPerChunk) !=0)
247 0 : chunksize(nchunk-1)=wConvSize%planesPerChunk;
248 :
249 0 : Int warner=0;
250 0 : Vector<Complex> maxes(wConvSize);
251 : // Bool maxdel;
252 : // Complex* maxptr=maxes.getStorage(maxdel);
253 0 : Matrix<Complex> corr(inner, inner);
254 0 : Vector<Complex> correction(inner);
255 0 : for (Int iy=-inner/2;iy<inner/2;iy++) {
256 :
257 0 : ggridder.correctX1D(correction, iy+inner/2);
258 0 : corr.row(iy+inner/2)=correction;
259 : }
260 : Bool cpcor;
261 :
262 0 : Complex *cor=corr.getStorage(cpcor);
263 0 : Vector<Int>pcsupp;
264 0 : pcsupp=convSupport;
265 : Bool delsupstor;
266 0 : Int* suppstor=pcsupp.getStorage(delsupstor);
267 0 : Double s1=sampling(1);
268 0 : Double s0=sampling(0);
269 : ///////////Por FFTPack
270 0 : Vector<Float> wsave(2*convSize*convSize+15);
271 0 : Int lsav=2*convSize*convSize+15;
272 : Bool wsavesave;
273 0 : Float *wsaveptr=wsave.getStorage(wsavesave);
274 : Int ier;
275 0 : FFTPack::cfft2i(convSize, convSize, wsaveptr, lsav, ier);
276 : //////////
277 0 : Matrix<Complex> screen(convSize, convSize);
278 0 : makeGWplane(screen, 0, s0, s1, wsaveptr, lsav, inner, cor, wScale);
279 0 : Float maxconv=abs(screen(0,0));
280 0 : for (Int chunkId=nchunk-1; chunkId >= 0; --chunkId){
281 : //cerr << "chunkId " << chunkId << endl;
282 0 : Cube<Complex> convFuncSect( maxConvSize/2-1, maxConvSize/2-1, chunksize(chunkId));
283 0 : convFuncSect.set(0.0);
284 0 : Bool convFuncStor=false;
285 0 : Complex *convFuncPtr=convFuncSect.getStorage(convFuncStor);
286 : //////openmp like to share reference param ...but i don't like to share
287 0 : Int cpConvSize=maxConvSize;
288 : //cerr << "orig convsize " << convSize << endl;
289 : // Int cpWConvSize=wConvSize;
290 0 : Double cpWscale=wScale;
291 0 : Int wstart=planesPerChunk*chunkId;
292 0 : Int wend=wstart+chunksize(chunkId)-1;
293 : #ifdef _OPENMP
294 0 : omp_set_nested(0);
295 : #endif
296 0 : #pragma omp parallel for default(none) firstprivate(/*cpWConvSize,*/ cpConvSize, convFuncPtr, s0, s1, wsaveptr, lsav, cor, inner, cpWscale, wstart, wend)
297 : for (Int iw=wstart; iw < (wend+1) ; ++iw) {
298 : Matrix<Complex> screen1(cpConvSize, cpConvSize);
299 : makeGWplane(screen1, iw, s0, s1, wsaveptr, lsav, inner, cor, cpWscale);
300 : ooLong offset=ooLong(ooLong(iw-wstart)*ooLong(cpConvSize/2-1)*ooLong(cpConvSize/2-1));
301 : for (ooLong y=0; y< ooLong(cpConvSize/2)-1; ++y){
302 : for (ooLong x=0; x< ooLong(cpConvSize/2)-1; ++x){
303 : ////////convFuncPtr[offset+y*(convSize/2-1)+x] = quarter(x,y);
304 : convFuncPtr[offset+y*ooLong(cpConvSize/2-1)+x] = screen1(x,y);
305 : }
306 : }
307 :
308 : }//iw
309 0 : convFuncSect.putStorage(convFuncPtr, convFuncStor);
310 0 : for (uInt iw=0; iw< uInt(chunksize(chunkId)); ++iw){
311 0 : convFuncSect.xyPlane(iw)=convFuncSect.xyPlane(iw)/(maxconv);
312 : }
313 0 : convFuncPtr=convFuncSect.getStorage(convFuncStor);
314 0 : Int thischunkSize=chunksize(chunkId);
315 : //cerr << "chunkId* planesPerChunk " << chunkId* planesPerChunk << " chunksize " << thischunkSize << endl;
316 0 : #pragma omp parallel for default(none) firstprivate(suppstor, cpConvSize, /*cpWConvSize,*/ convFuncPtr, maxConvSize, thischunkSize, wstart, planesPerChunk, chunkId) reduction(+: warner)
317 : for (Int iw=0; iw<thischunkSize; iw++) {
318 : Bool found=false;
319 : Int trial=0;
320 : ooLong ploffset=(ooLong)(cpConvSize/2-1)*(ooLong)(cpConvSize/2-1)*(ooLong)iw;
321 : for (trial=cpConvSize/2-2;trial>0;trial--) {
322 : // if((abs(convFunc(trial,0,iw))>1e-3)||(abs(convFunc(0,trial,iw))>1e-3) ) {
323 : if((abs(convFuncPtr[(ooLong)(trial)+ploffset])>1e-3)||(abs(convFuncPtr[(ooLong)(trial*(cpConvSize/2-1))+ploffset])>1e-3) ) {
324 : //cout <<"iw " << iw << " x " << abs(convFunc(trial,0,iw)) << " y "
325 : // <<abs(convFunc(0,trial,iw)) << endl;
326 : found=true;
327 : break;
328 : }
329 : }
330 : Int trueIw=iw+chunkId* planesPerChunk;
331 : if(found) {
332 : suppstor[trueIw]=Int(0.5+Float(trial)/Float(convSampling_p))+1;
333 : if(suppstor[trueIw]*convSampling_p*2 >= maxConvSize){
334 : suppstor[trueIw]=cpConvSize/2/convSampling_p-1;
335 : ++warner;
336 : }
337 : }
338 : else{
339 : suppstor[trueIw]=cpConvSize/2/convSampling_p-1;
340 : ++warner;
341 : }
342 : }
343 0 : convFuncSect.putStorage(convFuncPtr, convFuncStor);
344 0 : if(chunkId==(nchunk-1)){
345 0 : Int newConvSize=2*(suppstor[wConvSize-1]+2)*convSampling;
346 0 : if(newConvSize < convSize){
347 0 : convFunc.resize((newConvSize/2-1),
348 0 : (newConvSize/2-1),
349 : wConvSize);
350 0 : convSize=newConvSize;
351 : }
352 : else{
353 0 : convFunc.resize((convSize/2-1),
354 0 : (convSize/2-1),
355 : wConvSize);
356 : }
357 : }
358 0 : IPosition blc(3, 0,0,wstart);
359 0 : IPosition trc(3, (convSize/2-2),(convSize/2-2),wend);
360 0 : convFunc(blc, trc)=convFuncSect(IPosition(3, 0,0,0), IPosition(3, (convSize/2-2),
361 0 : (convSize/2-2), thischunkSize-1));
362 :
363 :
364 : }//chunkId
365 0 : corr.putStorage(cor, cpcor);
366 0 : pcsupp.putStorage(suppstor, delsupstor);
367 0 : convSupport=pcsupp;
368 0 : Double pbSum=0.0;
369 0 : for (Int iy=-convSupport(0);iy<=convSupport(0);iy++) {
370 0 : for (Int ix=-convSupport(0);ix<=convSupport(0);ix++) {
371 0 : pbSum+=real(convFunc(abs(ix)*convSampling_p,abs(iy)*convSampling_p,0));
372 : }
373 : }
374 0 : if(pbSum>0.0) {
375 0 : convFunc*=Complex(1.0/pbSum,0.0);
376 : }
377 : else {
378 : os << "Convolution function integral is not positive"
379 0 : << LogIO::EXCEPTION;
380 : }
381 :
382 :
383 0 : os << "Convolution support = " << convSupport*convSampling_p
384 : << " pixels in Fourier plane"
385 0 : << LogIO::POST;
386 0 : convSupportBlock_p.resize(actualConvIndex_p+1);
387 0 : convSupportBlock_p[actualConvIndex_p]= new Vector<Int>();
388 0 : convSupportBlock_p[actualConvIndex_p]->assign(convSupport);
389 0 : convFunctions_p.resize(actualConvIndex_p+1);
390 0 : convFunctions_p[actualConvIndex_p]= new Cube<Complex>(convFunc);
391 0 : convSizes_p.resize(actualConvIndex_p+1, true);
392 0 : convSizes_p(actualConvIndex_p)=convSize;
393 0 : Int maxMemoryMB=HostInfo::memoryTotal(true)/1024;
394 : Int memoryMB;
395 0 : memoryMB = Int(Double(convSize/2-1)*Double(convSize/2-1)*
396 0 : Double(wConvSize)*8.0/1024.0/1024.0);
397 : os << "Memory used in gridding function = "
398 : << memoryMB << " MB from maximum "
399 0 : << maxMemoryMB << " MB" << LogIO::POST;
400 0 : convSampling=convSampling_p;
401 0 : wScale=Float((wConvSize-1)*(wConvSize-1))/wScaler_p;
402 :
403 : }//end func
404 :
405 0 : void WPConvFunc::makeGWplane(Matrix<Complex>& screen, const Int iw, Double s0, Double s1, Float *& wsaveptr, Int& lsav, Int& inner, Complex*& cor, Double&cpWscale){
406 :
407 0 : Int cpConvSize=screen.shape()(0);
408 : //cerr << "cpConvSize " << cpConvSize << " shape " << screen.shape() << endl;
409 0 : screen.set(0.0);
410 : Bool cpscr;
411 0 : Complex *scr=screen.getStorage(cpscr);
412 0 : Double twoPiW=2.0*C::pi*Double(iw*iw)/cpWscale;
413 0 : for (Int iy=-inner/2;iy<inner/2;iy++) {
414 0 : Double m=s1*Double(iy);
415 0 : Double msq=m*m;
416 : //////Int offset= (iy+convSize/2)*convSize;
417 : ///fftpack likes it flipped
418 0 : ooLong offset= (iy>-1 ? ooLong(iy) : ooLong(iy+cpConvSize))*ooLong(cpConvSize);
419 0 : for (Int ix=-inner/2;ix<inner/2;ix++) {
420 : ////// Int ind=offset+ix+convSize/2;
421 : ///fftpack likes it flipped
422 0 : ooLong ind=offset+(ix > -1 ? ooLong(ix) : ooLong(ix+cpConvSize));
423 0 : Double l=s0*Double(ix);
424 0 : Double rsq=l*l+msq;
425 0 : if(rsq<1.0) {
426 0 : Double phase=twoPiW*(sqrt(1.0-rsq)-1.0);
427 : Double cval, sval;
428 0 : SINCOS(phase, sval, cval);
429 0 : Complex comval(cval, sval);
430 0 : scr[ind]=(cor[ix+inner/2+ (iy+inner/2)*inner])*comval;
431 : //screen(ix+convSize/2, iy+convSize/2)=comval;
432 : }
433 : }
434 : }
435 : ////Por FFTPack
436 0 : Vector<Float>work(2*cpConvSize*cpConvSize);
437 0 : Int lenwrk=2*cpConvSize*cpConvSize;
438 : Bool worksave;
439 0 : Float *workptr=work.getStorage(worksave);
440 : Int ier;
441 0 : FFTPack::cfft2f(cpConvSize, cpConvSize, cpConvSize, scr, wsaveptr, lsav, workptr, lenwrk, ier);
442 :
443 0 : screen.putStorage(scr, cpscr);
444 :
445 0 : }
446 0 : void WPConvFunc::findConvFunction2(const ImageInterface<Complex>& image,
447 : const VisBuffer& vb,
448 : const Int& wConvSizeUser,
449 : const Vector<Double>& uvScale,
450 : const Vector<Double>& uvOffset,
451 : const Float& padding,
452 : Int& convSampling,
453 : Cube<Complex>& convFunc,
454 : Int& convSize,
455 : Vector<Int>& convSupport, Double& wScale){
456 :
457 :
458 :
459 :
460 0 : if(checkCenterPix(image)){
461 0 : convFunc.resize();
462 0 : convFunc.reference(convFunc_p);
463 0 : convSize=convSize_p;
464 0 : convSampling=convSampling_p;
465 0 : convSupport.resize();
466 0 : convSupport=convSupport_p;
467 0 : wScale=Float((convFunc.shape()(2)-1)*(convFunc.shape()(2)-1))/wScaler_p;
468 0 : return;
469 : }
470 :
471 :
472 0 : LogIO os;
473 0 : os << LogOrigin("WPConvFunc", "findConvFunction") << LogIO::NORMAL;
474 :
475 :
476 : // Get the coordinate system
477 0 : CoordinateSystem coords(image.coordinates());
478 0 : Int directionIndex=coords.findCoordinate(Coordinate::DIRECTION);
479 0 : nx_p=Int(image.shape()(directionIndex));
480 0 : ny_p=Int(image.shape()(directionIndex+1));
481 :
482 0 : Int wConvSize=wConvSizeUser;
483 : ////Automatic mode
484 : Double maxUVW;
485 0 : if(wConvSize < 1){
486 : //cerr << "max, min, rms " << maxW_p << " " << minW_p << " " << rmsW_p << endl;
487 0 : maxUVW=rmsW_p < 0.5*(minW_p+maxW_p) ? 1.05*maxW_p: (rmsW_p /(0.5*((minW_p)+maxW_p))*1.05*maxW_p) ;
488 : //maxUVW=min(maxUVW, 0.25/abs(image.coordinates().increment()(0)));
489 0 : wConvSize=Int(maxUVW*fabs(sin(fabs(image.coordinates().increment()(0))*max(nx_p, ny_p)/2.0)));
490 : //maxUVW=1.05*maxW_p;
491 : //wConvSize=100*(maxUVW*maxUVW*image.coordinates().increment()(0)*image.coordinates().increment()(0));
492 : //cerr << "wConvSize 0 " << wConvSize << " nx_p " << nx_p << endl;
493 : //if(rmsW_p < 0.5*(minW_p+maxW_p))
494 : // wConvSize=wConvSize*(minW_p+maxW_p)*(minW_p+maxW_p)/(rmsW_p*rmsW_p);
495 :
496 0 : convSupport.resize(wConvSize);
497 : }
498 : else{
499 :
500 :
501 0 : if(maxW_p> 0.0)
502 0 : maxUVW= 1.05*maxW_p;
503 : else
504 0 : maxUVW=0.25/abs(image.coordinates().increment()(0));
505 :
506 : }
507 0 : if(wConvSize>1) {
508 0 : os << "W projection using " << wConvSize << " planes" << LogIO::POST;
509 :
510 : os << "Using maximum possible W = " << maxUVW
511 0 : << " (wavelengths)" << LogIO::POST;
512 :
513 0 : Double invLambdaC=vb.frequency()(0)/C::c;
514 : os << "Typical wavelength = " << 1.0/invLambdaC
515 0 : << " (m)" << LogIO::POST;
516 :
517 : // uvScale(2)=sqrt(Float(wConvSize-1))/maxUVW;
518 : // uvScale(2)=(Float(wConvSize-1))/maxUVW;
519 0 : wScale=Float((wConvSize-1)*(wConvSize-1))/maxUVW;
520 : //wScale=Float(wConvSize-1)/maxUVW;
521 0 : wScaler_p=maxUVW;;
522 0 : os << "Scaling in W (at maximum W) = " << 1.0/wScale
523 0 : << " wavelengths per pixel" << LogIO::POST;
524 : }
525 :
526 0 : Timer tim;
527 :
528 : // Set up the convolution function.
529 0 : if(wConvSize>1) {
530 : /* if(wConvSize>256) {
531 : convSampling=4;
532 : convSize=min(nx,ny);
533 : Int maxMemoryMB=HostInfo::memoryTotal()/1024;
534 : if(maxMemoryMB > 4000){
535 : convSize=min(convSize,1024);
536 : }
537 : else{
538 : convSize=min(convSize,512);
539 : }
540 :
541 : }
542 : else {
543 : convSampling=4;
544 : convSize=min(nx,ny);
545 : convSize=min(convSize,1024);
546 : }
547 : */
548 : // use memory size defined in aipsrc if exists
549 0 : Int maxMemoryMB=HostInfo::memoryTotal(true)/1024;
550 : //nominal 512 wprojplanes above that you may (or not) go swapping
551 0 : Double maxConvSizeConsidered=sqrt(Double(maxMemoryMB)/8.0*1024.0*1024.0/Double(wConvSize));
552 0 : CompositeNumber cn(Int(maxConvSizeConsidered/2.0)*2);
553 :
554 0 : convSampling_p=4;
555 0 : convSize=max(Int(nx_p*padding),Int(ny_p*padding));
556 0 : convSize=min(convSize,(Int)cn.nearestEven(Int(maxConvSizeConsidered/2.0)*2));
557 :
558 :
559 : }
560 : else {
561 0 : convSampling_p=1;
562 0 : convSize=max(Int(nx_p*padding),Int(ny_p*padding));
563 : }
564 0 : convSampling=convSampling_p;
565 0 : Int maxConvSize=convSize;
566 :
567 :
568 : // Make a two dimensional image to calculate the
569 : // primary beam. We want this on a fine grid in the
570 : // UV plane
571 0 : CompositeNumber cn(uInt(nx_p*2));
572 0 : AlwaysAssert(directionIndex>=0, AipsError);
573 0 : DirectionCoordinate dc=coords.directionCoordinate(directionIndex);
574 0 : Vector<Double> sampling;
575 0 : sampling = dc.increment();
576 0 : sampling*=Double(convSampling_p);
577 : //sampling*=Double(max(nx,ny))/Double(convSize);
578 0 : sampling[0]*=Double(cn.nextLargerEven(Int(padding*Float(nx_p)-0.5)))/Double(convSize);
579 0 : sampling[1]*=Double(cn.nextLargerEven(Int(padding*Float(ny_p)-0.5)))/Double(convSize);
580 0 : dc.setIncrement(sampling);
581 :
582 0 : Vector<Double> unitVec(2);
583 0 : unitVec=convSize/2;
584 0 : dc.setReferencePixel(unitVec);
585 :
586 : // Set the reference value to that of the image center for sure.
587 : {
588 : // dc.setReferenceValue(mTangent_p.getAngle().getValue());
589 0 : MDirection wcenter;
590 0 : Vector<Double> pcenter(2);
591 0 : pcenter(0) = nx_p/2;
592 0 : pcenter(1) = ny_p/2;
593 0 : coords.directionCoordinate(directionIndex).toWorld( wcenter, pcenter );
594 0 : dc.setReferenceValue(wcenter.getAngle().getValue());
595 : }
596 0 : coords.replaceCoordinate(dc, directionIndex);
597 : // coords.list(os, MDoppler::RADIO, IPosition(), IPosition());
598 :
599 0 : IPosition pbShape(4, convSize, convSize, 1, 1);
600 0 : TempImage<Complex> twoDPB(pbShape, coords);
601 :
602 0 : Int inner=convSize/convSampling_p;
603 : ConvolveGridder<Double, Complex>
604 0 : ggridder(IPosition(2, inner, inner), uvScale, uvOffset, "SF");
605 : /*
606 : ConvolveGridder<Double, Complex>
607 : ggridder(IPosition(2, cn.nextLargerEven(Int(padding*Float(nx_p)-0.5)),
608 : cn.nextLargerEven(Int(padding*Float(ny_p)-0.5))), uvScale,
609 : uvOffset, "SF");
610 : */
611 0 : convFunc.resize(); // break any reference
612 0 : convFunc.resize(convSize/2-1, convSize/2-1, wConvSize);
613 0 : convFunc.set(0.0);
614 0 : Bool convFuncStor=false;
615 0 : Complex *convFuncPtr=convFunc.getStorage(convFuncStor);
616 :
617 0 : IPosition start(4, 0, 0, 0, 0);
618 0 : IPosition pbSlice(4, convSize, convSize, 1, 1);
619 :
620 : //Bool writeResults=false;
621 0 : Int warner=0;
622 :
623 :
624 : // Accumulate terms
625 : //Matrix<Complex> screen(convSize, convSize);
626 : //screen.set(Complex(0.0));
627 0 : Vector<Complex> maxes(wConvSize);
628 : Bool maxdel;
629 0 : Complex* maxptr=maxes.getStorage(maxdel);
630 0 : Matrix<Complex> corr(inner, inner);
631 0 : Vector<Complex> correction(inner);
632 0 : for (Int iy=-inner/2;iy<inner/2;iy++) {
633 :
634 0 : ggridder.correctX1D(correction, iy+inner/2);
635 0 : corr.row(iy+inner/2)=correction;
636 : }
637 : Bool cpcor;
638 :
639 0 : Complex *cor=corr.getStorage(cpcor);
640 0 : Double s1=sampling(1);
641 0 : Double s0=sampling(0);
642 : ///////////Por FFTPack
643 0 : Vector<Float> wsave(2*convSize*convSize+15);
644 0 : Int lsav=2*convSize*convSize+15;
645 : Bool wsavesave;
646 0 : Float *wsaveptr=wsave.getStorage(wsavesave);
647 : Int ier;
648 0 : FFTPack::cfft2i(convSize, convSize, wsaveptr, lsav, ier);
649 : //////////
650 : #ifdef _OPENMP
651 0 : omp_set_nested(0);
652 : #endif
653 : //////openmp like to share reference param ...but i don't like to share
654 0 : Int cpConvSize=convSize;
655 0 : Int cpWConvSize=wConvSize;
656 0 : Double cpWscale=wScale;
657 : //Float max0=1.0;
658 0 : #pragma omp parallel for default(none) firstprivate(cpWConvSize, cpConvSize, convFuncPtr, s0, s1, wsaveptr, ier, lsav, cor, inner, maxptr, cpWscale, C::pi )
659 :
660 : for (Int iw=0; iw< cpWConvSize;iw++) {
661 : // First the w term
662 : Matrix<Complex> screen(cpConvSize, cpConvSize);
663 : screen=0.0;
664 : Bool cpscr;
665 : Complex *scr=screen.getStorage(cpscr);
666 : if(cpWConvSize>1) {
667 : // Double twoPiW=2.0*C::pi*sqrt(Double(iw))/uvScale(2);
668 : //Double twoPiW=2.0*C::pi*Double(iw)/wScale_p;
669 : Double twoPiW=2.0*C::pi*Double(iw*iw)/cpWscale;
670 : for (Int iy=-inner/2;iy<inner/2;iy++) {
671 : Double m=s1*Double(iy);
672 : Double msq=m*m;
673 : //////Int offset= (iy+convSize/2)*convSize;
674 : ///fftpack likes it flipped
675 : ooLong offset= (iy>-1 ? ooLong(iy) : ooLong(iy+cpConvSize))*ooLong(cpConvSize);
676 : for (Int ix=-inner/2;ix<inner/2;ix++) {
677 : ////// Int ind=offset+ix+convSize/2;
678 : ///fftpack likes it flipped
679 : ooLong ind=offset+(ix > -1 ? ooLong(ix) : ooLong(ix+cpConvSize));
680 : Double l=s0*Double(ix);
681 : Double rsq=l*l+msq;
682 : if(rsq<1.0) {
683 : Double phase=twoPiW*(sqrt(1.0-rsq)-1.0);
684 : Double cval, sval;
685 : SINCOS(phase, sval, cval);
686 : Complex comval(cval, sval);
687 : scr[ind]=(cor[ix+inner/2+ (iy+inner/2)*inner])*comval;
688 : //screen(ix+convSize/2, iy+convSize/2)=comval;
689 : }
690 : }
691 : }
692 : }
693 : else {
694 : screen=1.0;
695 : }
696 : /////////screen.putStorage(scr, cpscr);
697 :
698 : // spheroidal function
699 : /* Vector<Complex> correction(inner);
700 : for (Int iy=-inner/2;iy<inner/2;iy++) {
701 : ggridder.correctX1D(correction, iy+inner/2);
702 : for (Int ix=-inner/2;ix<inner/2;ix++) {
703 : screen(ix+convSize/2,iy+convSize/2)*=correction(ix+inner/2);
704 : }
705 : }
706 : */
707 : //twoDPB.putSlice(screen, IPosition(4, 0));
708 : // Write out screen as an image
709 : /*if(writeResults) {
710 : ostringstream name;
711 : name << "Screen" << iw+1;
712 : if(Table::canDeleteTable(name)) Table::deleteTable(name);
713 : PagedImage<Float> thisScreen(pbShape, coords, name);
714 : LatticeExpr<Float> le(real(twoDPB));
715 : thisScreen.copyData(le);
716 : }
717 : */
718 :
719 :
720 : // Now FFT and get the result back
721 : //LatticeFFT::cfft2d(twoDPB);
722 : /////////Por FFTPack
723 : Vector<Float>work(2*cpConvSize*cpConvSize);
724 : Int lenwrk=2*cpConvSize*cpConvSize;
725 : Bool worksave;
726 : Float *workptr=work.getStorage(worksave);
727 : FFTPack::cfft2f(cpConvSize, cpConvSize, cpConvSize, scr, wsaveptr, lsav, workptr, lenwrk, ier);
728 :
729 : screen.putStorage(scr, cpscr);
730 : /////////////////////
731 : // Write out FT of screen as an image
732 : /*if(1) {
733 : CoordinateSystem ftCoords(coords);
734 : directionIndex=ftCoords.findCoordinate(Coordinate::DIRECTION);
735 : AlwaysAssert(directionIndex>=0, AipsError);
736 : dc=coords.directionCoordinate(directionIndex);
737 : Vector<Bool> axes(2); axes(0)=true;axes(1)=true;
738 : Vector<Int> shape(2); shape(0)=convSize;shape(1)=convSize;
739 : Coordinate* ftdc=dc.makeFourierCoordinate(axes,shape);
740 : ftCoords.replaceCoordinate(*ftdc, directionIndex);
741 : delete ftdc; ftdc=0;
742 : ostringstream name;
743 : name << "FTScreen" << iw+1;
744 : if(Table::canDeleteTable(name)) Table::deleteTable(name);
745 : PagedImage<Float> thisScreen(pbShape, ftCoords, name);
746 : //LatticeExpr<Float> le(real(twoDPB));
747 : //thisScreen.copyData(le);
748 : thisScreen.put(real(screen));
749 : }
750 : */
751 : ////////IPosition start(4, convSize/2, convSize/2, 0, 0);
752 : ////////IPosition pbSlice(4, convSize/2-1, convSize/2-1, 1, 1);
753 : ///////convFunc.xyPlane(iw)=twoDPB.getSlice(start, pbSlice, true);
754 : //////Matrix<Complex> quarter(twoDPB.getSlice(start, pbSlice, true));
755 : // cerr << "quartershape " << quarter.shape() << endl;
756 : ooLong offset=ooLong(ooLong(iw)*ooLong(cpConvSize/2-1)*ooLong(cpConvSize/2-1));
757 : // cerr << "offset " << offset << " convfuncshape " << convFunc.shape() << " convSize " << convSize << endl;
758 : maxptr[iw]=screen(0,0);
759 : for (ooLong y=0; y< ooLong(cpConvSize/2)-1; ++y){
760 : for (ooLong x=0; x< ooLong(cpConvSize/2)-1; ++x){
761 : ////////convFuncPtr[offset+y*(convSize/2-1)+x] = quarter(x,y);
762 : convFuncPtr[offset+y*ooLong(cpConvSize/2-1)+x] = screen(x,y);
763 : }
764 : }
765 : }
766 :
767 0 : convFunc.putStorage(convFuncPtr, convFuncStor);
768 0 : corr.putStorage(cor, cpcor);
769 0 : maxes.putStorage(maxptr, maxdel);
770 0 : cerr << "maxes " << maxes << endl;
771 : //tim.show("After convFunc making ");
772 : //Complex maxconv=max(abs(convFunc));
773 0 : Complex maxconv=max(abs(maxes));
774 : //cerr << maxes << " maxconv " << maxconv << endl;
775 : //Do it by plane as the / operator makes a copy of the whole array
776 0 : for (uInt iw=0; iw< uInt(wConvSize); ++iw)
777 0 : convFunc.xyPlane(iw)=convFunc.xyPlane(iw)/real(maxconv);
778 : //tim.show("After convFunc norming ");
779 : // Find the edge of the function by stepping in from the
780 : // uv plane edge. We do this for each plane to save time on the
781 : // gridding (about a factor of two)
782 0 : convSupport=-1;
783 0 : Vector<Int> pcsupp=convSupport;
784 : #ifdef _OPENMP
785 0 : omp_set_nested(0);
786 : #endif
787 : Bool delsupstor;
788 0 : Int* suppstor=pcsupp.getStorage(delsupstor);
789 0 : convFuncPtr=convFunc.getStorage(convFuncStor);
790 0 : #pragma omp parallel for default(none) firstprivate(suppstor, cpConvSize, cpWConvSize, convFuncPtr, maxConvSize) reduction(+: warner)
791 : for (Int iw=0;iw<cpWConvSize;iw++) {
792 : Bool found=false;
793 : Int trial=0;
794 : ooLong ploffset=(ooLong)(cpConvSize/2-1)*(ooLong)(cpConvSize/2-1)*(ooLong)iw;
795 : for (trial=cpConvSize/2-2;trial>0;trial--) {
796 : // if((abs(convFunc(trial,0,iw))>1e-3)||(abs(convFunc(0,trial,iw))>1e-3) ) {
797 : if((abs(convFuncPtr[(ooLong)(trial)+ploffset])>1e-3)||(abs(convFuncPtr[(ooLong)(trial*(cpConvSize/2-1))+ploffset])>1e-3) ) {
798 : //cout <<"iw " << iw << " x " << abs(convFunc(trial,0,iw)) << " y "
799 : // <<abs(convFunc(0,trial,iw)) << endl;
800 : found=true;
801 : break;
802 : }
803 : }
804 : if(found) {
805 : suppstor[iw]=Int(0.5+Float(trial)/Float(convSampling_p))+1;
806 : if(suppstor[iw]*convSampling_p*2 >= maxConvSize){
807 : suppstor[iw]=cpConvSize/2/convSampling_p-1;
808 : ++warner;
809 : }
810 : }
811 : else{
812 : suppstor[iw]=cpConvSize/2/convSampling_p-1;
813 : ++warner;
814 : }
815 : }
816 0 : pcsupp.putStorage(suppstor, delsupstor);
817 0 : convSupport=pcsupp;
818 : //tim.show("After suppport locing ");
819 0 : if(convSupport(0)<1) {
820 : os << "Convolution function is misbehaved - support seems to be zero"
821 0 : << LogIO::EXCEPTION;
822 : }
823 :
824 0 : if(warner > 5) {
825 : os << LogIO::WARN
826 : <<"Many of the Convolution functions go beyond " << maxConvSize
827 0 : <<" pixels allocated" << LogIO::POST;
828 : os << LogIO::WARN
829 : << "You may consider reducing the size of your image or use facets"
830 0 : << LogIO::POST;
831 : }
832 : /*
833 : if(1) {
834 : CoordinateSystem ftCoords(coords);
835 : Int directionIndex=ftCoords.findCoordinate(Coordinate::DIRECTION);
836 : AlwaysAssert(directionIndex>=0, AipsError);
837 : dc=coords.directionCoordinate(directionIndex);
838 : Vector<Bool> axes(2); axes(0)=true;axes(1)=true;
839 : Vector<Int> shape(2); shape(0)=convSize;shape(1)=convSize;
840 : Coordinate* ftdc=dc.makeFourierCoordinate(axes,shape);
841 : ftCoords.replaceCoordinate(*ftdc, directionIndex);
842 : delete ftdc; ftdc=0;
843 : ostringstream name;
844 : name << "FTScreenWproj" ;
845 : if(Table::canDeleteTable(name)) Table::deleteTable(name);
846 : PagedImage<Complex> thisScreen(IPosition(4, convFunc.shape()(0), convFunc.shape()(1), 1, convFunc.shape()(2)), ftCoords, name);
847 : thisScreen.put(convFunc.reform(IPosition(4, convFunc.shape()(0), convFunc.shape()(1), 1, convFunc.shape()(2))));
848 : }
849 : */
850 :
851 :
852 : // Normalize such that plane 0 sums to 1 (when jumping in
853 : // steps of convSampling)
854 0 : Double pbSum=0.0;
855 0 : for (Int iy=-convSupport(0);iy<=convSupport(0);iy++) {
856 0 : for (Int ix=-convSupport(0);ix<=convSupport(0);ix++) {
857 0 : pbSum+=real(convFunc(abs(ix)*convSampling_p,abs(iy)*convSampling_p,0));
858 : }
859 : }
860 0 : if(pbSum>0.0) {
861 0 : convFunc*=Complex(1.0/pbSum,0.0);
862 : }
863 : else {
864 : os << "Convolution function integral is not positive"
865 0 : << LogIO::EXCEPTION;
866 : }
867 0 : os << "Convolution support = " << convSupport*convSampling_p
868 : << " pixels in Fourier plane"
869 0 : << LogIO::POST;
870 :
871 : //tim.show("After pbsumming ");
872 :
873 0 : convSupportBlock_p.resize(actualConvIndex_p+1);
874 0 : convSupportBlock_p[actualConvIndex_p]= new Vector<Int>();
875 0 : convSupportBlock_p[actualConvIndex_p]->assign(convSupport);
876 0 : convFunctions_p.resize(actualConvIndex_p+1);
877 0 : convFunctions_p[actualConvIndex_p]= new Cube<Complex>();
878 0 : Int newConvSize=2*(max(convSupport)+2)*convSampling;
879 :
880 0 : if(newConvSize < convSize){
881 0 : IPosition blc(3, 0,0,0);
882 0 : IPosition trc(3, (newConvSize/2-2),
883 0 : (newConvSize/2-2),
884 0 : convSupport.shape()(0)-1);
885 :
886 0 : *(convFunctions_p[actualConvIndex_p])=convFunc(blc,trc);
887 : // convFunctions_p[actualConvIndex_p]->assign(Cube<Complex>(convFunc(blc,trc)));
888 0 : convSize=newConvSize;
889 : }
890 : else{
891 0 : *(convFunctions_p[actualConvIndex_p])=convFunc;
892 : }
893 : // read out memory size from aisprc if exists
894 0 : Int maxMemoryMB=HostInfo::memoryTotal(true)/1024;
895 : Int memoryMB;
896 0 : memoryMB = Int(Double(convSize/2-1)*Double(convSize/2-1)*
897 0 : Double(wConvSize)*8.0/1024.0/1024.0);
898 : os << "Memory used in gridding function = "
899 : << memoryMB << " MB from maximum "
900 0 : << maxMemoryMB << " MB" << LogIO::POST;
901 0 : convFunc.resize();
902 0 : convFunc.reference(*convFunctions_p[actualConvIndex_p]);
903 0 : convSizes_p.resize(actualConvIndex_p+1, true);
904 0 : convSizes_p(actualConvIndex_p)=convSize;
905 :
906 0 : convSampling=convSampling_p;
907 0 : wScale=Float((wConvSize-1)*(wConvSize-1))/wScaler_p;
908 : //tim.show("After calculating WConv funx ");
909 :
910 :
911 :
912 : }
913 :
914 0 : Bool WPConvFunc::checkCenterPix(const ImageInterface<Complex>& image){
915 :
916 0 : CoordinateSystem imageCoord=image.coordinates();
917 0 : MDirection wcenter;
918 0 : Int directionIndex=imageCoord.findCoordinate(Coordinate::DIRECTION);
919 : DirectionCoordinate
920 0 : directionCoord=imageCoord.directionCoordinate(directionIndex);
921 0 : Vector<Double> incr=directionCoord.increment();
922 0 : nx_p=image.shape()(directionIndex);
923 0 : ny_p=image.shape()(directionIndex+1);
924 :
925 :
926 : //Images with same number of pixels and increments can have the same conv functions
927 0 : ostringstream oos;
928 0 : oos << std::setprecision(6);
929 :
930 0 : oos << nx_p << "_"<< fabs(incr(0)) << "_";
931 0 : oos << ny_p << "_"<< fabs(incr(1));
932 0 : String imageKey(oos);
933 :
934 0 : if(convFunctionMap_p.size( ) == 0){
935 0 : convFunctionMap_p.insert(std::pair<casacore::String, casacore::Int>(imageKey, 0));
936 0 : actualConvIndex_p=0;
937 0 : return false;
938 : }
939 :
940 0 : if( convFunctionMap_p.find(imageKey) == convFunctionMap_p.end( ) ){
941 0 : actualConvIndex_p=convFunctionMap_p.size( );
942 0 : convFunctionMap_p.insert(std::pair<casacore::String, casacore::Int>(imageKey,actualConvIndex_p));
943 0 : return false;
944 : }
945 : else{
946 0 : actualConvIndex_p=convFunctionMap_p[imageKey];
947 0 : convFunc_p.resize(); // break any reference
948 0 : convFunc_p.reference(*convFunctions_p[actualConvIndex_p]);
949 0 : convSupport_p.resize();
950 0 : convSupport_p.reference(*convSupportBlock_p[actualConvIndex_p]);
951 0 : convSize_p=convSizes_p[actualConvIndex_p];
952 :
953 : }
954 :
955 0 : return true;
956 : }
957 :
958 0 : Bool WPConvFunc::toRecord(RecordInterface& rec){
959 :
960 0 : Int numConv=convFunctions_p.nelements();
961 : try{
962 0 : rec.define("numconv", numConv);
963 0 : auto convptr = convFunctionMap_p.begin( );
964 0 : for (Int k=0; k < numConv; ++k, ++convptr){
965 0 : rec.define("convfunctions"+String::toString(k), *(convFunctions_p[k]));
966 0 : rec.define("convsupportblock"+String::toString(k), *(convSupportBlock_p[k]));
967 0 : rec.define("key"+String::toString(k), convptr->first);
968 0 : rec.define("val"+String::toString(k), convptr->second);
969 : }
970 0 : rec.define("convsizes", convSizes_p);
971 0 : rec.define("actualconvIndex",actualConvIndex_p);
972 0 : rec.define("convsize", convSize_p);
973 0 : rec.define("convsupport", convSupport_p);
974 0 : rec.define("convfunc",convFunc_p);
975 0 : rec.define("wscaler", wScaler_p);
976 0 : rec.define("convsampling", convSampling_p);
977 0 : rec.define("nx", nx_p);
978 0 : rec.define("ny", ny_p);
979 : }
980 0 : catch(AipsError x) {
981 0 : return false;
982 : }
983 0 : return true;
984 :
985 :
986 :
987 : }
988 :
989 0 : Bool WPConvFunc::fromRecord(String& err, const RecordInterface& rec){
990 :
991 0 : Int numConv=0;
992 : try{
993 0 : rec.get("numconv", numConv);
994 0 : convFunctions_p.resize(numConv, true, false);
995 0 : convSupportBlock_p.resize(numConv, true, false);
996 0 : convFunctionMap_p.clear( );
997 0 : for (Int k=0; k < numConv; ++k){
998 0 : convFunctions_p[k]=new Cube<Complex>();
999 0 : convSupportBlock_p[k]=new Vector<Int>();
1000 0 : rec.get("convfunctions"+String::toString(k), *(convFunctions_p[k]));
1001 0 : rec.get("convsupportblock"+String::toString(k), *(convSupportBlock_p[k]));
1002 0 : String key;
1003 : Int val;
1004 0 : rec.get("key"+String::toString(k), key);
1005 0 : rec.get("val"+String::toString(k), val);
1006 0 : convFunctionMap_p.insert(std::pair<casacore::String, casacore::Int>(key,val));
1007 : }
1008 0 : rec.get("convsizes", convSizes_p);
1009 0 : rec.get("actualconvIndex",actualConvIndex_p);
1010 0 : rec.get("convsize", convSize_p);
1011 0 : rec.get("convsupport", convSupport_p);
1012 0 : rec.get("convfunc",convFunc_p);
1013 0 : if(rec.isDefined("wscaler"))
1014 0 : rec.get("wscaler", wScaler_p);
1015 0 : rec.get("convsampling", convSampling_p);
1016 0 : rec.get("nx", nx_p);
1017 0 : rec.get("ny", ny_p);
1018 : }
1019 0 : catch(AipsError x) {
1020 0 : err=x.getMesg();
1021 0 : return false;
1022 : }
1023 0 : return true;
1024 :
1025 : }
1026 :
1027 :
1028 :
1029 :
1030 :
1031 : } //# NAMESPACE CASA - END
|
