Line data Source code
1 : //# WPConvFunc.cc: implementation of WPConvFunc
2 : //# Copyright (C) 2007-2016
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 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 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/Utilities/Assert.h>
40 : #include <casacore/casa/Utilities/CompositeNumber.h>
41 : #include <casacore/coordinates/Coordinates/CoordinateSystem.h>
42 : #include <casacore/coordinates/Coordinates/DirectionCoordinate.h>
43 :
44 : #include <casacore/images/Images/ImageInterface.h>
45 : #include <casacore/images/Images/PagedImage.h>
46 : #include <casacore/images/Images/TempImage.h>
47 : #include <casacore/casa/Logging/LogIO.h>
48 : #include <casacore/casa/Logging/LogSink.h>
49 : #include <casacore/casa/Logging/LogMessage.h>
50 :
51 : #include <casacore/lattices/Lattices/ArrayLattice.h>
52 : #include <casacore/lattices/Lattices/SubLattice.h>
53 : #include <casacore/lattices/LRegions/LCBox.h>
54 : #include <casacore/lattices/LEL/LatticeExpr.h>
55 : #include <casacore/lattices/Lattices/LatticeCache.h>
56 : #include <casacore/lattices/LatticeMath/LatticeFFT.h>
57 : #include <casacore/scimath/Mathematics/ConvolveGridder.h>
58 : #include <msvis/MSVis/VisBuffer2.h>
59 : #include <msvis/MSVis/VisibilityIterator2.h>
60 : #include <casacore/scimath/Mathematics/FFTPack.h>
61 : #include <msvis/MSVis/VisBuffer.h>
62 : #include <msvis/MSVis/VisibilityIterator.h>
63 : #include <synthesis/TransformMachines/SimplePBConvFunc.h> //por SINCOS
64 :
65 : #include <synthesis/TransformMachines2/WPConvFunc.h>
66 : #ifdef _OPENMP
67 : #include <omp.h>
68 : #endif
69 :
70 :
71 :
72 :
73 : using namespace casacore;
74 : namespace casa { //# NAMESPACE CASA - BEGIN
75 : namespace refim{ //namespace for refactoring imager
76 :
77 : using namespace casacore;
78 : using namespace casa;
79 : using namespace casacore;
80 : using namespace casa::refim;
81 : typedef unsigned long long ooLong;
82 0 : WPConvFunc::WPConvFunc(const Double minW, const Double maxW, const Double rmsW):
83 0 : actualConvIndex_p(-1), convSize_p(0), convSupport_p(0), minW_p(minW), maxW_p(maxW), rmsW_p(rmsW) {
84 : //
85 0 : }
86 0 : WPConvFunc::WPConvFunc(const WPConvFunc& other):
87 0 : actualConvIndex_p(-1), convSize_p(0), convSupport_p(0){
88 :
89 0 : operator=(other);
90 0 : }
91 :
92 0 : WPConvFunc& WPConvFunc::operator=(const WPConvFunc& other){
93 0 : if(this != &other){
94 0 : uInt numConv=other.convFunctions_p.nelements();
95 0 : convFunctions_p.resize(numConv, true, false);
96 0 : convSupportBlock_p.resize(numConv, true, false);
97 0 : for (uInt k=0; k < numConv; ++k){
98 0 : convFunctions_p[k]=new Cube<Complex>();
99 0 : *(convFunctions_p[k])=*(other.convFunctions_p[k]);
100 0 : convSupportBlock_p[k]=new Vector<Int> ();
101 0 : *(convSupportBlock_p[k]) = *(other.convSupportBlock_p[k]);
102 : }
103 :
104 0 : convFunctionMap_p=other.convFunctionMap_p;
105 0 : convSizes_p.resize();
106 0 : convSizes_p=other.convSizes_p;
107 0 : actualConvIndex_p=other.actualConvIndex_p;
108 0 : convSize_p=other.convSize_p;
109 0 : convSupport_p.resize();
110 0 : convSupport_p=other.convSupport_p;
111 0 : convFunc_p.resize();
112 0 : convFunc_p=other.convFunc_p;
113 0 : wScaler_p=other.wScaler_p;
114 0 : convSampling_p=other.convSampling_p;
115 0 : nx_p=other.nx_p;
116 0 : ny_p=other.ny_p;
117 0 : minW_p=other.minW_p;
118 0 : maxW_p=other.maxW_p;
119 0 : rmsW_p=other.rmsW_p;
120 :
121 :
122 :
123 : }
124 0 : return *this;
125 : }
126 :
127 0 : WPConvFunc::~WPConvFunc(){
128 : //usage of CountedPtr keeps this simple
129 :
130 0 : }
131 :
132 0 : WPConvFunc::WPConvFunc(const RecordInterface& rec):
133 0 : actualConvIndex_p(-1), convSize_p(0), convSupport_p(0){
134 0 : String error;
135 0 : if (!fromRecord(error, rec)) {
136 0 : throw (AipsError("Failed to create WPConvFunc: " + error));
137 : }
138 :
139 0 : }
140 :
141 0 : void WPConvFunc::findConvFunction(const ImageInterface<Complex>& image,
142 : const vi::VisBuffer2& vb,
143 : const Int& wConvSizeUser,
144 : const Vector<Double>& uvScale,
145 : const Vector<Double>& uvOffset,
146 : const Float& padding,
147 : Int& convSampling,
148 : Cube<Complex>& convFunc,
149 : Int& convSize,
150 : Vector<Int>& convSupport, Double& wScale){
151 0 : if(checkCenterPix(image)){
152 0 : convFunc.resize();
153 0 : convFunc.reference(convFunc_p);
154 0 : convSize=convSize_p;
155 0 : convSampling=convSampling_p;
156 0 : convSupport.resize();
157 0 : convSupport=convSupport_p;
158 0 : wScale=Float((convFunc.shape()(2)-1)*(convFunc.shape()(2)-1))/wScaler_p;
159 0 : return;
160 : }
161 0 : LogIO os;
162 0 : os << LogOrigin("WPConvFunc", "findConvFunction") << LogIO::NORMAL;
163 :
164 :
165 : // Get the coordinate system
166 0 : CoordinateSystem coords(image.coordinates());
167 0 : Int directionIndex=coords.findCoordinate(Coordinate::DIRECTION);
168 0 : nx_p=Int(image.shape()(directionIndex));
169 0 : ny_p=Int(image.shape()(directionIndex+1));
170 :
171 0 : Int wConvSize=wConvSizeUser;
172 : ////Automatic mode
173 : Double maxUVW;
174 0 : if(wConvSize < 1){
175 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) ;
176 0 : wConvSize=Int(maxUVW*fabs(sin(fabs(image.coordinates().increment()(0))*max(nx_p, ny_p)/2.0)));
177 0 : convSupport.resize(wConvSize);
178 : }
179 : else{
180 0 : if(maxW_p> 0.0)
181 0 : maxUVW= 1.05*maxW_p;
182 : else
183 0 : maxUVW=0.25/abs(image.coordinates().increment()(0));
184 :
185 : }
186 0 : if(wConvSize>1) {
187 0 : os << "W projection using " << wConvSize << " planes" << LogIO::POST;
188 :
189 : os << "Using maximum possible W = " << maxUVW
190 0 : << " (wavelengths)" << LogIO::POST;
191 :
192 0 : Double invLambdaC=vb.getFrequency(0,0)/C::c;
193 : os << "Typical wavelength = " << 1.0/invLambdaC
194 0 : << " (m)" << LogIO::POST;
195 :
196 :
197 0 : wScale=Float((wConvSize-1)*(wConvSize-1))/maxUVW;
198 : //wScale=Float(wConvSize-1)/maxUVW;
199 0 : wScaler_p=maxUVW;;
200 0 : os << "Scaling in W (at maximum W) = " << 1.0/wScale
201 0 : << " wavelengths per pixel" << LogIO::POST;
202 : }
203 :
204 0 : Int planesPerChunk=100;
205 :
206 :
207 0 : if(wConvSize>1) {
208 0 : Int maxMemoryMB=min(uInt64(HostInfo::memoryTotal(true)), uInt64(HostInfo::memoryFree()))/1024;
209 : ////Common you have to have 4 GB...memoryFree sometimes is whacky
210 0 : if(maxMemoryMB < 4000)
211 0 : maxMemoryMB=4000;
212 0 : convSize=max(Int(nx_p*padding),Int(ny_p*padding));
213 : ///Do the same thing as in WProject::init
214 0 : CompositeNumber cn(convSize);
215 0 : convSize = cn.nextLargerEven(convSize);
216 : //nominal 100 wprojplanes above that you may (or not) go swapping
217 :
218 0 : planesPerChunk=Int((Double(maxMemoryMB)/8.0*1024.0*1024.0)/Double(convSize)/Double(convSize));
219 : //cerr << "planesPerChunk" << planesPerChunk << endl;
220 0 : planesPerChunk=min(planesPerChunk, 100);
221 : // CompositeNumber cn(Int(maxConvSizeConsidered/2.0)*2);
222 0 : convSampling_p=4;
223 :
224 : //convSize=min(convSize,(Int)cn.nearestEven(Int(maxConvSizeConsidered/2.0)*2));
225 :
226 : }
227 : else{
228 0 : convSampling_p=1;
229 0 : convSize=max(Int(nx_p*padding),Int(ny_p*padding));
230 : }
231 0 : convSampling=convSampling_p;
232 0 : Int maxConvSize=convSize;
233 0 : convSupport.resize(wConvSize);
234 0 : convSupport.set(-1);
235 : ////set sampling
236 0 : AlwaysAssert(directionIndex>=0, AipsError);
237 0 : DirectionCoordinate dc=coords.directionCoordinate(directionIndex);
238 0 : Vector<Double> sampling;
239 0 : sampling = dc.increment();
240 0 : sampling*=Double(convSampling_p);
241 : //sampling*=Double(max(nx,ny))/Double(convSize);
242 0 : sampling[0]*=Double(padding*Float(nx_p))/Double(convSize);
243 0 : sampling[1]*=Double(padding*Float(ny_p))/Double(convSize);
244 : /////
245 0 : Int inner=convSize/convSampling_p;
246 : ConvolveGridder<Double, Complex>
247 0 : ggridder(IPosition(2, inner, inner), uvScale, uvOffset, "SF");
248 :
249 0 : Int nchunk=wConvSize/planesPerChunk;
250 0 : if((wConvSize%planesPerChunk) !=0)
251 0 : nchunk +=1;
252 0 : Vector<Int> chunksize(nchunk, planesPerChunk);
253 0 : if((wConvSize%planesPerChunk) !=0)
254 0 : chunksize(nchunk-1)=wConvSize%planesPerChunk;
255 :
256 0 : Int warner=0;
257 0 : Vector<Complex> maxes(wConvSize);
258 : // Bool maxdel;
259 : // Complex* maxptr=maxes.getStorage(maxdel);
260 0 : Matrix<Complex> corr(inner, inner);
261 0 : Vector<Complex> correction(inner);
262 0 : for (Int iy=-inner/2;iy<inner/2;iy++) {
263 :
264 0 : ggridder.correctX1D(correction, iy+inner/2);
265 0 : corr.row(iy+inner/2)=correction;
266 : }
267 : Bool cpcor;
268 :
269 0 : Complex *cor=corr.getStorage(cpcor);
270 0 : Vector<Int>pcsupp;
271 0 : pcsupp=convSupport;
272 : Bool delsupstor;
273 0 : Int* suppstor=pcsupp.getStorage(delsupstor);
274 0 : Double s1=sampling(1);
275 0 : Double s0=sampling(0);
276 : ///////////Por FFTPack
277 0 : Vector<Float> wsave(2*convSize*convSize+15);
278 0 : Int lsav=2*convSize*convSize+15;
279 : Bool wsavesave;
280 0 : Float *wsaveptr=wsave.getStorage(wsavesave);
281 : Int ier;
282 0 : FFTPack::cfft2i(convSize, convSize, wsaveptr, lsav, ier);
283 : //////////
284 0 : Matrix<Complex> screen(convSize, convSize);
285 0 : makeGWplane(screen, 0, s0, s1, wsaveptr, lsav, inner, cor, wScale);
286 0 : Float maxconv=abs(screen(0,0));
287 0 : for (Int chunkId=nchunk-1; chunkId >= 0; --chunkId){
288 : //cerr << "chunkId " << chunkId << endl;
289 0 : Cube<Complex> convFuncSect( maxConvSize/2-1, maxConvSize/2-1, chunksize(chunkId));
290 0 : convFuncSect.set(0.0);
291 0 : Bool convFuncStor=false;
292 0 : Complex *convFuncPtr=convFuncSect.getStorage(convFuncStor);
293 : //////openmp like to share reference param ...but i don't like to share
294 0 : Int cpConvSize=maxConvSize;
295 : //cerr << "orig convsize " << convSize << endl;
296 : // Int cpWConvSize=wConvSize;
297 0 : Double cpWscale=wScale;
298 0 : Int wstart=planesPerChunk*chunkId;
299 0 : Int wend=wstart+chunksize(chunkId)-1;
300 :
301 : #ifdef _OPENMP
302 0 : omp_set_nested(0);
303 0 : Int nthreads=omp_get_max_threads();
304 0 : nthreads=min(nthreads, chunksize(chunkId));
305 0 : omp_set_num_threads(nthreads);
306 :
307 : #endif
308 0 : #pragma omp parallel for default(none) firstprivate(/*cpWConvSize,*/ cpConvSize, convFuncPtr, s0, s1, wsaveptr, lsav, cor, inner, cpWscale, wstart, wend) schedule(dynamic, 1)
309 :
310 :
311 : for (Int iw=wstart; iw < (wend+1) ; ++iw) {
312 : Matrix<Complex> screen1(cpConvSize, cpConvSize);
313 : makeGWplane(screen1, iw, s0, s1, wsaveptr, lsav, inner, cor, cpWscale);
314 : ooLong offset=ooLong(ooLong(iw-wstart)*ooLong(cpConvSize/2-1)*ooLong(cpConvSize/2-1));
315 : for (ooLong y=0; y< ooLong(cpConvSize/2)-1; ++y){
316 : for (ooLong x=0; x< ooLong(cpConvSize/2)-1; ++x){
317 : ////////convFuncPtr[offset+y*(convSize/2-1)+x] = quarter(x,y);
318 : convFuncPtr[offset+y*ooLong(cpConvSize/2-1)+x] = screen1(x,y);
319 : }
320 : }
321 :
322 : }//iw
323 :
324 0 : convFuncSect.putStorage(convFuncPtr, convFuncStor);
325 0 : for (uInt iw=0; iw< uInt(chunksize(chunkId)); ++iw){
326 0 : convFuncSect.xyPlane(iw)=convFuncSect.xyPlane(iw)/(maxconv);
327 : }
328 0 : convFuncPtr=convFuncSect.getStorage(convFuncStor);
329 0 : Int thischunkSize=chunksize(chunkId);
330 : //cerr << "chunkId* planesPerChunk " << chunkId* planesPerChunk << " chunksize " << thischunkSize << endl;
331 0 : #pragma omp parallel for default(none) firstprivate(suppstor, cpConvSize, /*cpWConvSize,*/ convFuncPtr, maxConvSize, thischunkSize, wstart, planesPerChunk, chunkId) reduction(+: warner)
332 : for (Int iw=0; iw<thischunkSize; iw++) {
333 : Bool found=false;
334 : Int trial=0;
335 : ooLong ploffset=(ooLong)(cpConvSize/2-1)*(ooLong)(cpConvSize/2-1)*(ooLong)iw;
336 : for (trial=cpConvSize/2-2;trial>0;trial--) {
337 : // if((abs(convFunc(trial,0,iw))>1e-3)||(abs(convFunc(0,trial,iw))>1e-3) ) {
338 : if((abs(convFuncPtr[(ooLong)(trial)+ploffset])>1e-3)||(abs(convFuncPtr[(ooLong)(trial*(cpConvSize/2-1))+ploffset])>1e-3) ) {
339 : //cout <<"iw " << iw << " x " << abs(convFunc(trial,0,iw)) << " y "
340 : // <<abs(convFunc(0,trial,iw)) << endl;
341 : found=true;
342 : break;
343 : }
344 : }
345 : Int trueIw=iw+chunkId* planesPerChunk;
346 : if(found) {
347 : suppstor[trueIw]=Int(0.5+Float(trial)/Float(convSampling_p))+1;
348 : if(suppstor[trueIw]*convSampling_p*2 >= maxConvSize){
349 : suppstor[trueIw]=cpConvSize/2/convSampling_p-1;
350 : ++warner;
351 : }
352 : }
353 : else{
354 : suppstor[trueIw]=cpConvSize/2/convSampling_p-1;
355 : ++warner;
356 : }
357 : }
358 0 : convFuncSect.putStorage(convFuncPtr, convFuncStor);
359 0 : if(chunkId==(nchunk-1)){
360 0 : Int newConvSize=2*(suppstor[wConvSize-1]+2)*convSampling;
361 0 : if(newConvSize < convSize){
362 0 : convFunc.resize((newConvSize/2-1),
363 0 : (newConvSize/2-1),
364 : wConvSize);
365 0 : convSize=newConvSize;
366 : }
367 : else{
368 0 : convFunc.resize((convSize/2-1),
369 0 : (convSize/2-1),
370 : wConvSize);
371 : }
372 : }
373 0 : IPosition blc(3, 0,0,wstart);
374 0 : IPosition trc(3, (convSize/2-2),(convSize/2-2),wend);
375 0 : convFunc(blc, trc)=convFuncSect(IPosition(3, 0,0,0), IPosition(3, (convSize/2-2),
376 0 : (convSize/2-2), thischunkSize-1));
377 :
378 :
379 : }//chunkId
380 0 : corr.putStorage(cor, cpcor);
381 0 : pcsupp.putStorage(suppstor, delsupstor);
382 0 : convSupport=pcsupp;
383 0 : Double pbSum=0.0;
384 0 : for (Int iy=-convSupport(0);iy<=convSupport(0);iy++) {
385 0 : for (Int ix=-convSupport(0);ix<=convSupport(0);ix++) {
386 0 : pbSum+=real(convFunc(abs(ix)*convSampling_p,abs(iy)*convSampling_p,0));
387 : }
388 : }
389 0 : if(pbSum>0.0) {
390 0 : convFunc*=Complex(1.0/pbSum,0.0);
391 : }
392 : else {
393 : os << "Convolution function integral is not positive"
394 0 : << LogIO::EXCEPTION;
395 : }
396 :
397 :
398 0 : os << "Convolution support = " << convSupport*convSampling_p
399 : << " pixels in Fourier plane"
400 0 : << LogIO::POST;
401 0 : convSupportBlock_p.resize(actualConvIndex_p+1);
402 0 : convSupportBlock_p[actualConvIndex_p]= new Vector<Int>();
403 0 : convSupportBlock_p[actualConvIndex_p]->assign(convSupport);
404 0 : convFunctions_p.resize(actualConvIndex_p+1);
405 0 : convFunctions_p[actualConvIndex_p]= new Cube<Complex>(convFunc);
406 0 : convSizes_p.resize(actualConvIndex_p+1, true);
407 0 : convSizes_p(actualConvIndex_p)=convSize;
408 0 : Int maxMemoryMB=HostInfo::memoryTotal(true)/1024;
409 : Int memoryMB;
410 0 : memoryMB = Int(Double(convSize/2-1)*Double(convSize/2-1)*
411 0 : Double(wConvSize)*8.0/1024.0/1024.0);
412 : os << "Memory used in gridding function = "
413 : << memoryMB << " MB from maximum "
414 0 : << maxMemoryMB << " MB" << LogIO::POST;
415 0 : convSampling=convSampling_p;
416 0 : wScale=Float((wConvSize-1)*(wConvSize-1))/wScaler_p;
417 :
418 : }//end func
419 :
420 0 : void WPConvFunc::makeGWplane(Matrix<Complex>& screen, const Int iw, Double s0, Double s1, Float *& wsaveptr, Int& lsav, Int& inner, Complex*& cor, Double&cpWscale){
421 :
422 0 : Int cpConvSize=screen.shape()(0);
423 : //cerr << "cpConvSize " << cpConvSize << " shape " << screen.shape() << endl;
424 0 : screen.set(0.0);
425 : Bool cpscr;
426 0 : Complex *scr=screen.getStorage(cpscr);
427 0 : Double twoPiW=2.0*C::pi*Double(iw*iw)/cpWscale;
428 0 : for (Int iy=-inner/2;iy<inner/2;iy++) {
429 0 : Double m=s1*Double(iy);
430 0 : Double msq=m*m;
431 : //////Int offset= (iy+convSize/2)*convSize;
432 : ///fftpack likes it flipped
433 0 : ooLong offset= (iy>-1 ? ooLong(iy) : ooLong(iy+cpConvSize))*ooLong(cpConvSize);
434 0 : for (Int ix=-inner/2;ix<inner/2;ix++) {
435 : ////// Int ind=offset+ix+convSize/2;
436 : ///fftpack likes it flipped
437 0 : ooLong ind=offset+(ix > -1 ? ooLong(ix) : ooLong(ix+cpConvSize));
438 0 : Double l=s0*Double(ix);
439 0 : Double rsq=l*l+msq;
440 0 : if(rsq<1.0) {
441 0 : Double phase=twoPiW*(sqrt(1.0-rsq)-1.0);
442 : Double cval, sval;
443 0 : SINCOS(phase, sval, cval);
444 0 : Complex comval(cval, sval);
445 0 : scr[ind]=(cor[ix+inner/2+ (iy+inner/2)*inner])*comval;
446 : //screen(ix+convSize/2, iy+convSize/2)=comval;
447 : }
448 : }
449 : }
450 : ////Por FFTPack
451 0 : Vector<Float>work(2*cpConvSize*cpConvSize);
452 0 : Int lenwrk=2*cpConvSize*cpConvSize;
453 : Bool worksave;
454 0 : Float *workptr=work.getStorage(worksave);
455 : Int ier;
456 0 : FFTPack::cfft2f(cpConvSize, cpConvSize, cpConvSize, scr, wsaveptr, lsav, workptr, lenwrk, ier);
457 :
458 0 : screen.putStorage(scr, cpscr);
459 :
460 0 : }
461 :
462 :
463 0 : Bool WPConvFunc::checkCenterPix(const ImageInterface<Complex>& image){
464 :
465 0 : CoordinateSystem imageCoord=image.coordinates();
466 0 : MDirection wcenter;
467 0 : Int directionIndex=imageCoord.findCoordinate(Coordinate::DIRECTION);
468 : DirectionCoordinate
469 0 : directionCoord=imageCoord.directionCoordinate(directionIndex);
470 0 : Vector<Double> incr=directionCoord.increment();
471 0 : nx_p=image.shape()(directionIndex);
472 0 : ny_p=image.shape()(directionIndex+1);
473 :
474 :
475 : //Images with same number of pixels and increments can have the same conv functions
476 0 : ostringstream oos;
477 0 : oos << std::setprecision(6);
478 :
479 0 : oos << nx_p << "_"<< fabs(incr(0)) << "_";
480 0 : oos << ny_p << "_"<< fabs(incr(1));
481 0 : String imageKey(oos);
482 :
483 0 : if(convFunctionMap_p.size( ) == 0){
484 0 : convFunctionMap_p.insert(std::pair<casacore::String, casacore::Int>(imageKey, 0));
485 0 : actualConvIndex_p=0;
486 0 : return false;
487 : }
488 :
489 0 : if( convFunctionMap_p.find(imageKey) == convFunctionMap_p.end( ) ){
490 0 : actualConvIndex_p=convFunctionMap_p.size( );
491 0 : convFunctionMap_p.insert(std::pair<casacore::String, casacore::Int>(imageKey,actualConvIndex_p));
492 0 : return false;
493 : }
494 : else{
495 0 : actualConvIndex_p=convFunctionMap_p[imageKey];
496 0 : convFunc_p.resize(); // break any reference
497 0 : convFunc_p.reference(*convFunctions_p[actualConvIndex_p]);
498 0 : convSupport_p.resize();
499 0 : convSupport_p.reference(*convSupportBlock_p[actualConvIndex_p]);
500 0 : convSize_p=convSizes_p[actualConvIndex_p];
501 :
502 : }
503 :
504 0 : return true;
505 : }
506 :
507 0 : Bool WPConvFunc::toRecord(RecordInterface& rec){
508 :
509 0 : Int numConv=convFunctions_p.nelements();
510 : try{
511 0 : rec.define("numconv", numConv);
512 0 : auto convptr = convFunctionMap_p.begin( );
513 0 : for (Int k=0; k < numConv; ++k, ++convptr){
514 0 : rec.define("convfunctions"+String::toString(k), *(convFunctions_p[k]));
515 0 : rec.define("convsupportblock"+String::toString(k), *(convSupportBlock_p[k]));
516 0 : rec.define("key"+String::toString(k), convptr->first);
517 0 : rec.define("val"+String::toString(k), convptr->second);
518 : }
519 0 : rec.define("convsizes", convSizes_p);
520 0 : rec.define("actualconvIndex",actualConvIndex_p);
521 0 : rec.define("convsize", convSize_p);
522 0 : rec.define("convsupport", convSupport_p);
523 0 : rec.define("convfunc",convFunc_p);
524 0 : rec.define("wscaler", wScaler_p);
525 0 : rec.define("convsampling", convSampling_p);
526 0 : rec.define("nx", nx_p);
527 0 : rec.define("ny", ny_p);
528 : }
529 0 : catch(AipsError x) {
530 0 : return false;
531 : }
532 0 : return true;
533 :
534 :
535 :
536 : }
537 :
538 0 : Bool WPConvFunc::fromRecord(String& err, const RecordInterface& rec){
539 :
540 0 : Int numConv=0;
541 : try{
542 0 : rec.get("numconv", numConv);
543 0 : convFunctions_p.resize(numConv, true, false);
544 0 : convSupportBlock_p.resize(numConv, true, false);
545 0 : convFunctionMap_p.clear( );
546 0 : for (Int k=0; k < numConv; ++k){
547 0 : convFunctions_p[k]=new Cube<Complex>();
548 0 : convSupportBlock_p[k]=new Vector<Int>();
549 0 : rec.get("convfunctions"+String::toString(k), *(convFunctions_p[k]));
550 0 : rec.get("convsupportblock"+String::toString(k), *(convSupportBlock_p[k]));
551 0 : String key;
552 : Int val;
553 0 : rec.get("key"+String::toString(k), key);
554 0 : rec.get("val"+String::toString(k), val);
555 0 : convFunctionMap_p.insert(std::pair<casacore::String, casacore::Int>(key,val));
556 : }
557 0 : rec.get("convsizes", convSizes_p);
558 0 : rec.get("actualconvIndex",actualConvIndex_p);
559 0 : rec.get("convsize", convSize_p);
560 0 : rec.get("convsupport", convSupport_p);
561 0 : rec.get("convfunc",convFunc_p);
562 0 : if(rec.isDefined("wscaler"))
563 0 : rec.get("wscaler", wScaler_p);
564 0 : rec.get("convsampling", convSampling_p);
565 0 : rec.get("nx", nx_p);
566 0 : rec.get("ny", ny_p);
567 : }
568 0 : catch(AipsError x) {
569 0 : err=x.getMesg();
570 0 : return false;
571 : }
572 0 : return true;
573 :
574 : }
575 :
576 :
577 :
578 :
579 : } //end of namespace refim
580 : } //# NAMESPACE CASA - END
|
