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: casa-feedback@nrao.edu.
21 : //# Postal address: AIPS++ Project Office
22 : //# National Radio Astronomy Observatory
23 : //# 520 Edgemont Road
24 : //# Charlottesville, VA 22903-2475 USA
25 : //#
26 : //#
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 91 : WPConvFunc::WPConvFunc(const Double minW, const Double maxW, const Double rmsW):
83 91 : actualConvIndex_p(-1), convSize_p(0), convSupport_p(0), minW_p(minW), maxW_p(maxW), rmsW_p(rmsW) {
84 : //
85 91 : }
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 180 : WPConvFunc::~WPConvFunc(){
128 : //usage of CountedPtr keeps this simple
129 :
130 180 : }
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 1824 : 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 1824 : if(checkCenterPix(image)){
152 1787 : convFunc.resize();
153 1787 : convFunc.reference(convFunc_p);
154 1787 : convSize=convSize_p;
155 1787 : convSampling=convSampling_p;
156 1787 : convSupport.resize();
157 1787 : convSupport=convSupport_p;
158 1787 : wScale=Float((convFunc.shape()(2)-1)*(convFunc.shape()(2)-1))/wScaler_p;
159 1787 : return;
160 : }
161 37 : LogIO os;
162 37 : os << LogOrigin("WPConvFunc", "findConvFunction") << LogIO::NORMAL;
163 :
164 :
165 : // Get the coordinate system
166 37 : CoordinateSystem coords(image.coordinates());
167 37 : Int directionIndex=coords.findCoordinate(Coordinate::DIRECTION);
168 37 : nx_p=Int(image.shape()(directionIndex));
169 37 : ny_p=Int(image.shape()(directionIndex+1));
170 :
171 37 : Int wConvSize=wConvSizeUser;
172 : ////Automatic mode
173 : Double maxUVW;
174 37 : 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 32 : if(maxW_p> 0.0)
181 0 : maxUVW= 1.05*maxW_p;
182 : else
183 32 : maxUVW=0.25/abs(image.coordinates().increment()(0));
184 :
185 : }
186 37 : if(wConvSize>1) {
187 37 : os << "W projection using " << wConvSize << " planes" << LogIO::POST;
188 :
189 : os << "Using maximum possible W = " << maxUVW
190 37 : << " (wavelengths)" << LogIO::POST;
191 :
192 37 : Double invLambdaC=vb.getFrequency(0,0)/C::c;
193 : os << "Typical wavelength = " << 1.0/invLambdaC
194 37 : << " (m)" << LogIO::POST;
195 :
196 :
197 37 : wScale=Float((wConvSize-1)*(wConvSize-1))/maxUVW;
198 : //wScale=Float(wConvSize-1)/maxUVW;
199 37 : wScaler_p=maxUVW;;
200 37 : os << "Scaling in W (at maximum W) = " << 1.0/wScale
201 37 : << " wavelengths per pixel" << LogIO::POST;
202 : }
203 :
204 37 : Int planesPerChunk=100;
205 :
206 :
207 37 : if(wConvSize>1) {
208 37 : Int maxMemoryMB=min(uInt64(HostInfo::memoryTotal(true)), uInt64(HostInfo::memoryFree()))/1024;
209 : ////Common you have to have 4 GB...memoryFree sometimes is whacky
210 37 : if(maxMemoryMB < 4000)
211 0 : maxMemoryMB=4000;
212 37 : convSize=max(Int(nx_p*padding),Int(ny_p*padding));
213 : ///Do the same thing as in WProject::init
214 37 : CompositeNumber cn(convSize);
215 37 : convSize = cn.nextLargerEven(convSize);
216 : //nominal 100 wprojplanes above that you may (or not) go swapping
217 :
218 37 : planesPerChunk=Int((Double(maxMemoryMB)/8.0*1024.0*1024.0)/Double(convSize)/Double(convSize));
219 : //cerr << "planesPerChunk" << planesPerChunk << endl;
220 37 : planesPerChunk=min(planesPerChunk, 100);
221 : // CompositeNumber cn(Int(maxConvSizeConsidered/2.0)*2);
222 37 : convSampling_p=4;
223 :
224 : //convSize=min(convSize,(Int)cn.nearestEven(Int(maxConvSizeConsidered/2.0)*2));
225 :
226 37 : }
227 : else{
228 0 : convSampling_p=1;
229 0 : convSize=max(Int(nx_p*padding),Int(ny_p*padding));
230 : }
231 37 : convSampling=convSampling_p;
232 37 : Int maxConvSize=convSize;
233 37 : convSupport.resize(wConvSize);
234 37 : convSupport.set(-1);
235 : ////set sampling
236 37 : AlwaysAssert(directionIndex>=0, AipsError);
237 37 : DirectionCoordinate dc=coords.directionCoordinate(directionIndex);
238 37 : Vector<Double> sampling;
239 37 : sampling = dc.increment();
240 37 : sampling*=Double(convSampling_p);
241 : //sampling*=Double(max(nx,ny))/Double(convSize);
242 37 : sampling[0]*=Double(padding*Float(nx_p))/Double(convSize);
243 37 : sampling[1]*=Double(padding*Float(ny_p))/Double(convSize);
244 : /////
245 37 : Int inner=convSize/convSampling_p;
246 : ConvolveGridder<Double, Complex>
247 74 : ggridder(IPosition(2, inner, inner), uvScale, uvOffset, "SF");
248 :
249 37 : Int nchunk=wConvSize/planesPerChunk;
250 37 : if((wConvSize%planesPerChunk) !=0)
251 37 : nchunk +=1;
252 37 : Vector<Int> chunksize(nchunk, planesPerChunk);
253 37 : if((wConvSize%planesPerChunk) !=0)
254 37 : chunksize(nchunk-1)=wConvSize%planesPerChunk;
255 :
256 37 : Int warner=0;
257 37 : Vector<Complex> maxes(wConvSize);
258 : // Bool maxdel;
259 : // Complex* maxptr=maxes.getStorage(maxdel);
260 37 : Matrix<Complex> corr(inner, inner);
261 37 : Vector<Complex> correction(inner);
262 10335 : for (Int iy=-inner/2;iy<inner/2;iy++) {
263 :
264 10298 : ggridder.correctX1D(correction, iy+inner/2);
265 10298 : corr.row(iy+inner/2)=correction;
266 : }
267 : Bool cpcor;
268 :
269 37 : Complex *cor=corr.getStorage(cpcor);
270 37 : Vector<Int>pcsupp;
271 37 : pcsupp=convSupport;
272 : Bool delsupstor;
273 37 : Int* suppstor=pcsupp.getStorage(delsupstor);
274 37 : Double s1=sampling(1);
275 37 : Double s0=sampling(0);
276 : ///////////Por FFTPack
277 37 : Vector<Float> wsave(2*convSize*convSize+15);
278 37 : Int lsav=2*convSize*convSize+15;
279 : Bool wsavesave;
280 37 : Float *wsaveptr=wsave.getStorage(wsavesave);
281 : Int ier;
282 37 : FFTPack::cfft2i(convSize, convSize, wsaveptr, lsav, ier);
283 : //////////
284 37 : Matrix<Complex> screen(convSize, convSize);
285 37 : makeGWplane(screen, 0, s0, s1, wsaveptr, lsav, inner, cor, wScale);
286 37 : Float maxconv=abs(screen(0,0));
287 83 : for (Int chunkId=nchunk-1; chunkId >= 0; --chunkId){
288 : //cerr << "chunkId " << chunkId << endl;
289 46 : Cube<Complex> convFuncSect( maxConvSize/2-1, maxConvSize/2-1, chunksize(chunkId));
290 46 : convFuncSect.set(0.0);
291 46 : Bool convFuncStor=false;
292 46 : Complex *convFuncPtr=convFuncSect.getStorage(convFuncStor);
293 : //////openmp like to share reference param ...but i don't like to share
294 46 : Int cpConvSize=maxConvSize;
295 : //cerr << "orig convsize " << convSize << endl;
296 : // Int cpWConvSize=wConvSize;
297 46 : Double cpWscale=wScale;
298 46 : Int wstart=planesPerChunk*chunkId;
299 46 : Int wend=wstart+chunksize(chunkId)-1;
300 :
301 : #ifdef _OPENMP
302 46 : omp_set_nested(0);
303 46 : Int nthreads=omp_get_max_threads();
304 46 : nthreads=min(nthreads, chunksize(chunkId));
305 46 : omp_set_num_threads(nthreads);
306 :
307 : #endif
308 46 : #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 46 : convFuncSect.putStorage(convFuncPtr, convFuncStor);
325 1403 : for (uInt iw=0; iw< uInt(chunksize(chunkId)); ++iw){
326 1357 : convFuncSect.xyPlane(iw)=convFuncSect.xyPlane(iw)/(maxconv);
327 : }
328 46 : convFuncPtr=convFuncSect.getStorage(convFuncStor);
329 46 : Int thischunkSize=chunksize(chunkId);
330 : //cerr << "chunkId* planesPerChunk " << chunkId* planesPerChunk << " chunksize " << thischunkSize << endl;
331 46 : #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 46 : convFuncSect.putStorage(convFuncPtr, convFuncStor);
359 46 : if(chunkId==(nchunk-1)){
360 37 : Int newConvSize=2*(suppstor[wConvSize-1]+2)*convSampling;
361 37 : if(newConvSize < convSize){
362 37 : convFunc.resize((newConvSize/2-1),
363 37 : (newConvSize/2-1),
364 : wConvSize);
365 37 : convSize=newConvSize;
366 : }
367 : else{
368 0 : convFunc.resize((convSize/2-1),
369 0 : (convSize/2-1),
370 : wConvSize);
371 : }
372 : }
373 46 : IPosition blc(3, 0,0,wstart);
374 46 : IPosition trc(3, (convSize/2-2),(convSize/2-2),wend);
375 92 : convFunc(blc, trc)=convFuncSect(IPosition(3, 0,0,0), IPosition(3, (convSize/2-2),
376 92 : (convSize/2-2), thischunkSize-1));
377 :
378 :
379 46 : }//chunkId
380 37 : corr.putStorage(cor, cpcor);
381 37 : pcsupp.putStorage(suppstor, delsupstor);
382 37 : convSupport=pcsupp;
383 37 : Double pbSum=0.0;
384 370 : for (Int iy=-convSupport(0);iy<=convSupport(0);iy++) {
385 3330 : for (Int ix=-convSupport(0);ix<=convSupport(0);ix++) {
386 2997 : pbSum+=real(convFunc(abs(ix)*convSampling_p,abs(iy)*convSampling_p,0));
387 : }
388 : }
389 37 : if(pbSum>0.0) {
390 37 : 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 37 : os << "Convolution support = " << convSupport*convSampling_p
399 : << " pixels in Fourier plane"
400 37 : << LogIO::POST;
401 37 : convSupportBlock_p.resize(actualConvIndex_p+1);
402 37 : convSupportBlock_p[actualConvIndex_p]= new Vector<Int>();
403 37 : convSupportBlock_p[actualConvIndex_p]->assign(convSupport);
404 37 : convFunctions_p.resize(actualConvIndex_p+1);
405 37 : convFunctions_p[actualConvIndex_p]= new Cube<Complex>(convFunc);
406 37 : convSizes_p.resize(actualConvIndex_p+1, true);
407 37 : convSizes_p(actualConvIndex_p)=convSize;
408 37 : Int maxMemoryMB=HostInfo::memoryTotal(true)/1024;
409 : Int memoryMB;
410 37 : memoryMB = Int(Double(convSize/2-1)*Double(convSize/2-1)*
411 37 : Double(wConvSize)*8.0/1024.0/1024.0);
412 : os << "Memory used in gridding function = "
413 : << memoryMB << " MB from maximum "
414 37 : << maxMemoryMB << " MB" << LogIO::POST;
415 37 : convSampling=convSampling_p;
416 37 : wScale=Float((wConvSize-1)*(wConvSize-1))/wScaler_p;
417 :
418 37 : }//end func
419 :
420 1394 : 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 1394 : Int cpConvSize=screen.shape()(0);
423 : //cerr << "cpConvSize " << cpConvSize << " shape " << screen.shape() << endl;
424 1394 : screen.set(0.0);
425 : Bool cpscr;
426 1394 : Complex *scr=screen.getStorage(cpscr);
427 1394 : Double twoPiW=2.0*C::pi*Double(iw*iw)/cpWscale;
428 476672 : for (Int iy=-inner/2;iy<inner/2;iy++) {
429 475278 : Double m=s1*Double(iy);
430 475278 : Double msq=m*m;
431 : //////Int offset= (iy+convSize/2)*convSize;
432 : ///fftpack likes it flipped
433 475278 : ooLong offset= (iy>-1 ? ooLong(iy) : ooLong(iy+cpConvSize))*ooLong(cpConvSize);
434 234167482 : for (Int ix=-inner/2;ix<inner/2;ix++) {
435 : ////// Int ind=offset+ix+convSize/2;
436 : ///fftpack likes it flipped
437 233692204 : ooLong ind=offset+(ix > -1 ? ooLong(ix) : ooLong(ix+cpConvSize));
438 233692204 : Double l=s0*Double(ix);
439 233692204 : Double rsq=l*l+msq;
440 233692204 : if(rsq<1.0) {
441 233692204 : Double phase=twoPiW*(sqrt(1.0-rsq)-1.0);
442 : Double cval, sval;
443 233692204 : SINCOS(phase, sval, cval);
444 233692204 : Complex comval(cval, sval);
445 233692204 : 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 1394 : Vector<Float>work(2*cpConvSize*cpConvSize);
452 1394 : Int lenwrk=2*cpConvSize*cpConvSize;
453 : Bool worksave;
454 1394 : Float *workptr=work.getStorage(worksave);
455 : Int ier;
456 1394 : FFTPack::cfft2f(cpConvSize, cpConvSize, cpConvSize, scr, wsaveptr, lsav, workptr, lenwrk, ier);
457 :
458 1394 : screen.putStorage(scr, cpscr);
459 :
460 1394 : }
461 :
462 :
463 1824 : Bool WPConvFunc::checkCenterPix(const ImageInterface<Complex>& image){
464 :
465 1824 : CoordinateSystem imageCoord=image.coordinates();
466 1824 : MDirection wcenter;
467 1824 : Int directionIndex=imageCoord.findCoordinate(Coordinate::DIRECTION);
468 : DirectionCoordinate
469 1824 : directionCoord=imageCoord.directionCoordinate(directionIndex);
470 1824 : Vector<Double> incr=directionCoord.increment();
471 1824 : nx_p=image.shape()(directionIndex);
472 1824 : ny_p=image.shape()(directionIndex+1);
473 :
474 :
475 : //Images with same number of pixels and increments can have the same conv functions
476 1824 : ostringstream oos;
477 1824 : oos << std::setprecision(6);
478 :
479 1824 : oos << nx_p << "_"<< fabs(incr(0)) << "_";
480 1824 : oos << ny_p << "_"<< fabs(incr(1));
481 1824 : String imageKey(oos);
482 :
483 1824 : if(convFunctionMap_p.size( ) == 0){
484 37 : convFunctionMap_p.insert(std::pair<casacore::String, casacore::Int>(imageKey, 0));
485 37 : actualConvIndex_p=0;
486 37 : return false;
487 : }
488 :
489 1787 : 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 1787 : actualConvIndex_p=convFunctionMap_p[imageKey];
496 1787 : convFunc_p.resize(); // break any reference
497 1787 : convFunc_p.reference(*convFunctions_p[actualConvIndex_p]);
498 1787 : convSupport_p.resize();
499 1787 : convSupport_p.reference(*convSupportBlock_p[actualConvIndex_p]);
500 1787 : convSize_p=convSizes_p[actualConvIndex_p];
501 :
502 : }
503 :
504 1787 : return true;
505 1824 : }
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 0 : }
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 0 : }
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 0 : }
572 0 : return true;
573 :
574 : }
575 0 : Bool WPConvFunc::makeWConvFuncs(Cube<Complex>& wconv, Vector<Int>& supports, CoordinateSystem& cs, const Int& npix, const Vector<Double>& wVals){
576 0 : uInt nw=wVals.nelements();
577 0 : supports.resize(nw);
578 : //cerr << "wVals" << wVals << endl;
579 0 : Double wVal=wVals(nw-1);
580 0 : Matrix<Complex> screen;
581 0 : makeSkyWFunc(screen, cs, npix, wVal);
582 0 : FFT2D ft(True);
583 : Bool isCopyScr;
584 0 : Complex* scr=screen.getStorage(isCopyScr);
585 0 : ft.c2cFFT(scr, npix, npix);
586 0 : screen.putStorage(scr,isCopyScr);
587 0 : Int sup=findSupport(screen);
588 : ///TESTOO
589 : //sup=256;
590 : /////
591 0 : if (sup < 1)
592 0 : throw(AipsError("Could not find support for w term"));
593 0 : Int convSize=2*(sup+2);
594 0 : if (convSize > npix)
595 0 : convSize=npix;
596 0 : wconv.resize(convSize, convSize, nw);
597 0 : IPosition blc(2, npix/2-convSize/2, npix/2-convSize/2);
598 0 : IPosition trc(2, npix/2+convSize/2-1, npix/2+convSize/2-1);
599 0 : wconv.xyPlane(nw-1).assign(screen(blc, trc));
600 0 : supports(nw-1)=sup;
601 0 : for(Int k=0; k < (nw-1); ++k){
602 0 : wVal=wVals[k];
603 0 : makeSkyWFunc(screen, cs, npix, wVal);
604 0 : Complex* scr=screen.getStorage(isCopyScr);
605 0 : ft.c2cFFT(scr, npix, npix);
606 0 : screen.putStorage(scr,isCopyScr);
607 0 : sup=findSupport(screen);
608 0 : supports[k]=sup;
609 0 : wconv.xyPlane(k)=screen(blc, trc);
610 : }
611 :
612 :
613 :
614 0 : return True;
615 0 : }
616 0 : Int WPConvFunc::findSupport(Matrix<Complex>& arr){
617 0 : Int shape=arr.shape()[0];
618 0 : Float maxval=real(max((fabs(arr))));
619 0 : Float minval=real(min((fabs(arr))));
620 0 : Float maxshift=maxval-minval;
621 : //cerr << "####maxval " << maxval << " shift " << maxshift << endl;
622 0 : Int trial=0;
623 0 : Bool found=False;
624 0 : Float sqrt_2=1.0/sqrt(2);
625 0 : while(!found && (trial < (shape/2))){
626 0 : Float frac=(fabs( arr(shape/2, (shape/2+ trial) ) )-minval)/maxshift;
627 0 : if( frac < 1e-2){
628 0 : found=True;
629 : //cerr << "found at trial " << trial << endl;
630 :
631 : }
632 0 : ++trial;
633 :
634 : }
635 :
636 0 : if(!found)
637 0 : return shape/2;
638 :
639 0 : return trial;
640 :
641 : }
642 0 : Bool WPConvFunc::makeSkyWFunc(Matrix<Complex>& screen,
643 : const CoordinateSystem& cs, const Int& npix, const Double& wVal){
644 :
645 0 : screen.resize(npix, npix);
646 : //cerr << "SCR shape " << screen.shape() << endl;
647 0 : screen.set(0.0);
648 : //Assuming largest wVals is latest
649 0 : Double scale1=cs.increment()(0);
650 0 : Double scale2=cs.increment()(1);
651 0 : Double uvoff=Double(npix)/2.0;
652 : Bool isScrCopy;
653 0 : Complex *scr=screen.getStorage(isScrCopy);
654 :
655 0 : Double twoPiW=2.0*C::pi* wVal;
656 0 : for (uInt iy=0; iy < uInt(npix); ++iy){
657 0 : Double m=scale2*(Double(iy)-uvoff);
658 0 : Double msq=m*m;
659 0 : for(uInt ix=0; ix < uInt(npix); ++ix){
660 0 : Double l=scale1*(Double(ix)-uvoff);
661 0 : Double rsq=l*l+msq;
662 0 : if(rsq<1.0) {
663 0 : Double phase=twoPiW*(sqrt(1.0-rsq)-1.0);
664 : Double cval, sval;
665 0 : SINCOS(phase, sval, cval);
666 : //cerr << "ix " << ix << " iy " << iy << " pix " << ix+ iy*npix << endl;
667 0 : scr[ix+iy*npix]=Complex(cval, sval);
668 : //screen(ix, iy)=Complex(cval, sval);
669 : }
670 : }
671 :
672 :
673 : }
674 0 : screen.putStorage(scr, isScrCopy);
675 :
676 :
677 :
678 :
679 :
680 0 : return True;
681 :
682 :
683 :
684 : }
685 :
686 :
687 :
688 : } //end of namespace refim
689 : } //# NAMESPACE CASA - END
|
