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 80 : WPConvFunc::WPConvFunc(const Double minW, const Double maxW, const Double rmsW):
83 80 : actualConvIndex_p(-1), convSize_p(0), convSupport_p(0), minW_p(minW), maxW_p(maxW), rmsW_p(rmsW) {
84 : //
85 80 : }
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 160 : WPConvFunc::~WPConvFunc(){
128 : //usage of CountedPtr keeps this simple
129 :
130 160 : }
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 1544 : 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 1544 : if(checkCenterPix(image)){
152 1512 : convFunc.resize();
153 1512 : convFunc.reference(convFunc_p);
154 1512 : convSize=convSize_p;
155 1512 : convSampling=convSampling_p;
156 1512 : convSupport.resize();
157 1512 : convSupport=convSupport_p;
158 1512 : wScale=Float((convFunc.shape()(2)-1)*(convFunc.shape()(2)-1))/wScaler_p;
159 1512 : return;
160 : }
161 64 : LogIO os;
162 32 : os << LogOrigin("WPConvFunc", "findConvFunction") << LogIO::NORMAL;
163 :
164 :
165 : // Get the coordinate system
166 64 : CoordinateSystem coords(image.coordinates());
167 32 : Int directionIndex=coords.findCoordinate(Coordinate::DIRECTION);
168 32 : nx_p=Int(image.shape()(directionIndex));
169 32 : ny_p=Int(image.shape()(directionIndex+1));
170 :
171 32 : Int wConvSize=wConvSizeUser;
172 : ////Automatic mode
173 : Double maxUVW;
174 32 : if(wConvSize < 1){
175 5 : 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 5 : wConvSize=Int(maxUVW*fabs(sin(fabs(image.coordinates().increment()(0))*max(nx_p, ny_p)/2.0)));
177 5 : convSupport.resize(wConvSize);
178 : }
179 : else{
180 27 : if(maxW_p> 0.0)
181 0 : maxUVW= 1.05*maxW_p;
182 : else
183 27 : maxUVW=0.25/abs(image.coordinates().increment()(0));
184 :
185 : }
186 32 : if(wConvSize>1) {
187 32 : os << "W projection using " << wConvSize << " planes" << LogIO::POST;
188 :
189 : os << "Using maximum possible W = " << maxUVW
190 32 : << " (wavelengths)" << LogIO::POST;
191 :
192 32 : Double invLambdaC=vb.getFrequency(0,0)/C::c;
193 : os << "Typical wavelength = " << 1.0/invLambdaC
194 32 : << " (m)" << LogIO::POST;
195 :
196 :
197 32 : wScale=Float((wConvSize-1)*(wConvSize-1))/maxUVW;
198 : //wScale=Float(wConvSize-1)/maxUVW;
199 32 : wScaler_p=maxUVW;;
200 32 : os << "Scaling in W (at maximum W) = " << 1.0/wScale
201 32 : << " wavelengths per pixel" << LogIO::POST;
202 : }
203 :
204 32 : Int planesPerChunk=100;
205 :
206 :
207 32 : if(wConvSize>1) {
208 32 : Int maxMemoryMB=min(uInt64(HostInfo::memoryTotal(true)), uInt64(HostInfo::memoryFree()))/1024;
209 : ////Common you have to have 4 GB...memoryFree sometimes is whacky
210 32 : if(maxMemoryMB < 4000)
211 0 : maxMemoryMB=4000;
212 32 : convSize=max(Int(nx_p*padding),Int(ny_p*padding));
213 : ///Do the same thing as in WProject::init
214 32 : CompositeNumber cn(convSize);
215 32 : convSize = cn.nextLargerEven(convSize);
216 : //nominal 100 wprojplanes above that you may (or not) go swapping
217 :
218 32 : planesPerChunk=Int((Double(maxMemoryMB)/8.0*1024.0*1024.0)/Double(convSize)/Double(convSize));
219 : //cerr << "planesPerChunk" << planesPerChunk << endl;
220 32 : planesPerChunk=min(planesPerChunk, 100);
221 : // CompositeNumber cn(Int(maxConvSizeConsidered/2.0)*2);
222 32 : 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 32 : convSampling=convSampling_p;
232 32 : Int maxConvSize=convSize;
233 32 : convSupport.resize(wConvSize);
234 32 : convSupport.set(-1);
235 : ////set sampling
236 32 : AlwaysAssert(directionIndex>=0, AipsError);
237 64 : DirectionCoordinate dc=coords.directionCoordinate(directionIndex);
238 64 : Vector<Double> sampling;
239 32 : sampling = dc.increment();
240 32 : sampling*=Double(convSampling_p);
241 : //sampling*=Double(max(nx,ny))/Double(convSize);
242 32 : sampling[0]*=Double(padding*Float(nx_p))/Double(convSize);
243 32 : sampling[1]*=Double(padding*Float(ny_p))/Double(convSize);
244 : /////
245 32 : Int inner=convSize/convSampling_p;
246 : ConvolveGridder<Double, Complex>
247 96 : ggridder(IPosition(2, inner, inner), uvScale, uvOffset, "SF");
248 :
249 32 : Int nchunk=wConvSize/planesPerChunk;
250 32 : if((wConvSize%planesPerChunk) !=0)
251 32 : nchunk +=1;
252 64 : Vector<Int> chunksize(nchunk, planesPerChunk);
253 32 : if((wConvSize%planesPerChunk) !=0)
254 32 : chunksize(nchunk-1)=wConvSize%planesPerChunk;
255 :
256 32 : Int warner=0;
257 64 : Vector<Complex> maxes(wConvSize);
258 : // Bool maxdel;
259 : // Complex* maxptr=maxes.getStorage(maxdel);
260 64 : Matrix<Complex> corr(inner, inner);
261 64 : Vector<Complex> correction(inner);
262 9664 : for (Int iy=-inner/2;iy<inner/2;iy++) {
263 :
264 9632 : ggridder.correctX1D(correction, iy+inner/2);
265 9632 : corr.row(iy+inner/2)=correction;
266 : }
267 : Bool cpcor;
268 :
269 32 : Complex *cor=corr.getStorage(cpcor);
270 64 : Vector<Int>pcsupp;
271 32 : pcsupp=convSupport;
272 : Bool delsupstor;
273 32 : Int* suppstor=pcsupp.getStorage(delsupstor);
274 32 : Double s1=sampling(1);
275 32 : Double s0=sampling(0);
276 : ///////////Por FFTPack
277 64 : Vector<Float> wsave(2*convSize*convSize+15);
278 32 : Int lsav=2*convSize*convSize+15;
279 : Bool wsavesave;
280 32 : Float *wsaveptr=wsave.getStorage(wsavesave);
281 : Int ier;
282 32 : FFTPack::cfft2i(convSize, convSize, wsaveptr, lsav, ier);
283 : //////////
284 32 : Matrix<Complex> screen(convSize, convSize);
285 32 : makeGWplane(screen, 0, s0, s1, wsaveptr, lsav, inner, cor, wScale);
286 32 : Float maxconv=abs(screen(0,0));
287 73 : for (Int chunkId=nchunk-1; chunkId >= 0; --chunkId){
288 : //cerr << "chunkId " << chunkId << endl;
289 82 : Cube<Complex> convFuncSect( maxConvSize/2-1, maxConvSize/2-1, chunksize(chunkId));
290 41 : convFuncSect.set(0.0);
291 41 : Bool convFuncStor=false;
292 41 : Complex *convFuncPtr=convFuncSect.getStorage(convFuncStor);
293 : //////openmp like to share reference param ...but i don't like to share
294 41 : Int cpConvSize=maxConvSize;
295 : //cerr << "orig convsize " << convSize << endl;
296 : // Int cpWConvSize=wConvSize;
297 41 : Double cpWscale=wScale;
298 41 : Int wstart=planesPerChunk*chunkId;
299 41 : Int wend=wstart+chunksize(chunkId)-1;
300 :
301 : #ifdef _OPENMP
302 41 : omp_set_nested(0);
303 41 : Int nthreads=omp_get_max_threads();
304 41 : nthreads=min(nthreads, chunksize(chunkId));
305 41 : omp_set_num_threads(nthreads);
306 :
307 : #endif
308 41 : #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 41 : convFuncSect.putStorage(convFuncPtr, convFuncStor);
325 1332 : for (uInt iw=0; iw< uInt(chunksize(chunkId)); ++iw){
326 1291 : convFuncSect.xyPlane(iw)=convFuncSect.xyPlane(iw)/(maxconv);
327 : }
328 41 : convFuncPtr=convFuncSect.getStorage(convFuncStor);
329 41 : Int thischunkSize=chunksize(chunkId);
330 : //cerr << "chunkId* planesPerChunk " << chunkId* planesPerChunk << " chunksize " << thischunkSize << endl;
331 41 : #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 41 : convFuncSect.putStorage(convFuncPtr, convFuncStor);
359 41 : if(chunkId==(nchunk-1)){
360 32 : Int newConvSize=2*(suppstor[wConvSize-1]+2)*convSampling;
361 32 : if(newConvSize < convSize){
362 64 : convFunc.resize((newConvSize/2-1),
363 32 : (newConvSize/2-1),
364 : wConvSize);
365 32 : convSize=newConvSize;
366 : }
367 : else{
368 0 : convFunc.resize((convSize/2-1),
369 0 : (convSize/2-1),
370 : wConvSize);
371 : }
372 : }
373 82 : IPosition blc(3, 0,0,wstart);
374 41 : IPosition trc(3, (convSize/2-2),(convSize/2-2),wend);
375 82 : convFunc(blc, trc)=convFuncSect(IPosition(3, 0,0,0), IPosition(3, (convSize/2-2),
376 82 : (convSize/2-2), thischunkSize-1));
377 :
378 :
379 : }//chunkId
380 32 : corr.putStorage(cor, cpcor);
381 32 : pcsupp.putStorage(suppstor, delsupstor);
382 32 : convSupport=pcsupp;
383 32 : Double pbSum=0.0;
384 320 : for (Int iy=-convSupport(0);iy<=convSupport(0);iy++) {
385 2880 : for (Int ix=-convSupport(0);ix<=convSupport(0);ix++) {
386 2592 : pbSum+=real(convFunc(abs(ix)*convSampling_p,abs(iy)*convSampling_p,0));
387 : }
388 : }
389 32 : if(pbSum>0.0) {
390 32 : 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 32 : os << "Convolution support = " << convSupport*convSampling_p
399 : << " pixels in Fourier plane"
400 32 : << LogIO::POST;
401 32 : convSupportBlock_p.resize(actualConvIndex_p+1);
402 32 : convSupportBlock_p[actualConvIndex_p]= new Vector<Int>();
403 32 : convSupportBlock_p[actualConvIndex_p]->assign(convSupport);
404 32 : convFunctions_p.resize(actualConvIndex_p+1);
405 32 : convFunctions_p[actualConvIndex_p]= new Cube<Complex>(convFunc);
406 32 : convSizes_p.resize(actualConvIndex_p+1, true);
407 32 : convSizes_p(actualConvIndex_p)=convSize;
408 32 : Int maxMemoryMB=HostInfo::memoryTotal(true)/1024;
409 : Int memoryMB;
410 96 : memoryMB = Int(Double(convSize/2-1)*Double(convSize/2-1)*
411 64 : Double(wConvSize)*8.0/1024.0/1024.0);
412 : os << "Memory used in gridding function = "
413 : << memoryMB << " MB from maximum "
414 32 : << maxMemoryMB << " MB" << LogIO::POST;
415 32 : convSampling=convSampling_p;
416 32 : wScale=Float((wConvSize-1)*(wConvSize-1))/wScaler_p;
417 :
418 : }//end func
419 :
420 1323 : 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 1323 : Int cpConvSize=screen.shape()(0);
423 : //cerr << "cpConvSize " << cpConvSize << " shape " << screen.shape() << endl;
424 1323 : screen.set(0.0);
425 : Bool cpscr;
426 1323 : Complex *scr=screen.getStorage(cpscr);
427 1323 : Double twoPiW=2.0*C::pi*Double(iw*iw)/cpWscale;
428 472075 : for (Int iy=-inner/2;iy<inner/2;iy++) {
429 470752 : Double m=s1*Double(iy);
430 470752 : Double msq=m*m;
431 : //////Int offset= (iy+convSize/2)*convSize;
432 : ///fftpack likes it flipped
433 470752 : ooLong offset= (iy>-1 ? ooLong(iy) : ooLong(iy+cpConvSize))*ooLong(cpConvSize);
434 233616480 : for (Int ix=-inner/2;ix<inner/2;ix++) {
435 : ////// Int ind=offset+ix+convSize/2;
436 : ///fftpack likes it flipped
437 233145728 : ooLong ind=offset+(ix > -1 ? ooLong(ix) : ooLong(ix+cpConvSize));
438 233145728 : Double l=s0*Double(ix);
439 233145728 : Double rsq=l*l+msq;
440 233145728 : if(rsq<1.0) {
441 233145728 : Double phase=twoPiW*(sqrt(1.0-rsq)-1.0);
442 : Double cval, sval;
443 233145728 : SINCOS(phase, sval, cval);
444 233145728 : Complex comval(cval, sval);
445 233145728 : 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 2646 : Vector<Float>work(2*cpConvSize*cpConvSize);
452 1323 : Int lenwrk=2*cpConvSize*cpConvSize;
453 : Bool worksave;
454 1323 : Float *workptr=work.getStorage(worksave);
455 : Int ier;
456 1323 : FFTPack::cfft2f(cpConvSize, cpConvSize, cpConvSize, scr, wsaveptr, lsav, workptr, lenwrk, ier);
457 :
458 1323 : screen.putStorage(scr, cpscr);
459 :
460 1323 : }
461 :
462 :
463 1544 : Bool WPConvFunc::checkCenterPix(const ImageInterface<Complex>& image){
464 :
465 3088 : CoordinateSystem imageCoord=image.coordinates();
466 3088 : MDirection wcenter;
467 1544 : Int directionIndex=imageCoord.findCoordinate(Coordinate::DIRECTION);
468 : DirectionCoordinate
469 3088 : directionCoord=imageCoord.directionCoordinate(directionIndex);
470 3088 : Vector<Double> incr=directionCoord.increment();
471 1544 : nx_p=image.shape()(directionIndex);
472 1544 : ny_p=image.shape()(directionIndex+1);
473 :
474 :
475 : //Images with same number of pixels and increments can have the same conv functions
476 3088 : ostringstream oos;
477 1544 : oos << std::setprecision(6);
478 :
479 1544 : oos << nx_p << "_"<< fabs(incr(0)) << "_";
480 1544 : oos << ny_p << "_"<< fabs(incr(1));
481 3088 : String imageKey(oos);
482 :
483 1544 : if(convFunctionMap_p.size( ) == 0){
484 32 : convFunctionMap_p.insert(std::pair<casacore::String, casacore::Int>(imageKey, 0));
485 32 : actualConvIndex_p=0;
486 32 : return false;
487 : }
488 :
489 1512 : 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 1512 : actualConvIndex_p=convFunctionMap_p[imageKey];
496 1512 : convFunc_p.resize(); // break any reference
497 1512 : convFunc_p.reference(*convFunctions_p[actualConvIndex_p]);
498 1512 : convSupport_p.resize();
499 1512 : convSupport_p.reference(*convSupportBlock_p[actualConvIndex_p]);
500 1512 : convSize_p=convSizes_p[actualConvIndex_p];
501 :
502 : }
503 :
504 1512 : 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
|
