Line data Source code
1 : //# PSTerm.cc: implementation of PSTerm
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: 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 <synthesis/TransformMachines2/PSTerm.h>
29 : #include <synthesis/TransformMachines2/Utils.h>
30 : #include <synthesis/TransformMachines/SynthesisError.h>
31 : #ifdef _OPENMP
32 : #include <omp.h>
33 : #endif
34 :
35 : using namespace casacore;
36 : namespace casa { //# NAMESPACE CASA - BEGIN
37 : using namespace refim;
38 : //
39 : //-------------------------------------------------------------------------
40 : //
41 56 : void PSTerm::init(const IPosition shape,
42 : const Vector<Double>& uvScale,
43 : const Vector<Double>& uvOffset,
44 : const Double& psScale)
45 : {
46 56 : shape_p = shape;
47 56 : uvScale_p=uvScale;
48 56 : uvOffset_p = uvOffset;
49 56 : psScale_p = psScale;
50 :
51 56 : psCtor_p = new ConvolveGridder<Double, Complex>(shape_p, uvScale_p, uvOffset_p, "SF");
52 56 : sampling_p = (*psCtor_p).cSupport();
53 56 : support_p = (*psCtor_p).cSampling();
54 56 : }
55 0 : void PSTerm::reinit(const IPosition shape,
56 : const Vector<Double>& uvScale,
57 : const Vector<Double>& uvOffset,
58 : const Double& psScale)
59 : {
60 0 : if (shape.nelements() > 0) shape_p = shape;
61 0 : if (uvScale.nelements() > 0) uvScale_p=uvScale;
62 0 : if (uvOffset.nelements() > 0) uvOffset_p = uvOffset;
63 0 : if (psScale > 0) psScale_p = psScale;
64 :
65 0 : psCtor_p = new ConvolveGridder<Double, Complex>(shape_p, uvScale_p, uvOffset_p, "SF");
66 0 : sampling_p = (*psCtor_p).cSupport();
67 0 : support_p = (*psCtor_p).cSampling();
68 0 : }
69 : //
70 : //-------------------------------------------------------------------------
71 : //
72 0 : Matrix<Complex>& PSTerm::operator=(Matrix<Complex>& buf)
73 : {
74 0 : applySky(buf);
75 0 : return buf;
76 : }
77 : //
78 : //-------------------------------------------------------------------------
79 : //
80 0 : Matrix<Complex>& PSTerm::operator*=(Matrix<Complex>& buf)
81 : {
82 0 : applySky(buf,true);
83 0 : return buf;
84 : }
85 : //
86 : //-------------------------------------------------------------------------
87 : //
88 0 : void PSTerm::applySky(Matrix<Complex>& screen, Bool multiply)
89 : {
90 : // AlwaysAssert(psCtor_p.null()==false, SynthesisFTMachineError);
91 0 : Int nx=screen.shape()(0), ny=screen.shape()(1);
92 0 : Int convOrig=nx/2;
93 0 : Float xpart, psScale_local=psScale_p;
94 : #ifdef _OPENMP
95 0 : Int Nth=max(omp_get_max_threads()-2,1);
96 : #endif
97 :
98 0 : if (!isNoOp())
99 : {
100 0 : for (Int i=0; i<nx;i++)
101 : {
102 0 : xpart = square(i-convOrig);
103 : #ifdef _OPENMP
104 : //#pragma omp parallel default(none) firstprivate(xpart,convOrig, i) shared(screen,psScale_local,ny,multiply) num_threads(Nth)
105 0 : #pragma omp parallel firstprivate(xpart,convOrig, i) shared(psScale_local,ny,multiply) num_threads(Nth)
106 : #endif
107 : {
108 : #ifdef _OPENMP
109 : #pragma omp for
110 : #endif
111 : for (Int j=0;j<ny;j++)
112 : {
113 : Float ypart;
114 : ypart = sqrt(xpart + square(j-convOrig))*psScale_local;
115 : Float val = SynthesisUtils::libreSpheroidal(ypart)*(1.0-ypart*ypart);;
116 : if (val > 0)
117 : {
118 : if (multiply) screen(i, j) *= val;
119 : else screen(i, j) = val;
120 : }
121 : }
122 : }
123 : }
124 : }
125 0 : }
126 :
127 64 : void PSTerm::applySky(Matrix<Complex>& screen,
128 : const Vector<Double>& sampling,
129 : const Int inner)
130 : {
131 : //GG: convsize: 2048 inner: 512 uvscale: [-0.0397159, 0.0397159, 2.33547e-310] uvoffset:[1024, 1024, 0]
132 : // Vector<Double> uvScale(3);uvScale[0]=-0.0397159;uvScale[1]=0.0397159;uvScale[2]=0.0;
133 : // Vector<Double> uvOffset(3,0.0); uvOffset[0]=uvOffset[1]=1024;
134 : // psCtor_p=new ConvolveGridder<Double, Complex>(IPosition(2, inner, inner),
135 : // uvScale, uvOffset,
136 : // "SF");
137 :
138 :
139 64 : AlwaysAssert(psCtor_p.null()==false, SynthesisFTMachineError);
140 :
141 64 : Int convSize = screen.shape()[0];
142 : // Vector<Complex> correction(inner);
143 64 : Vector<Complex> correction;
144 64 : if (!isNoOp())
145 : {
146 32832 : for (Int iy=-inner/2;iy<inner/2;iy++)
147 : {
148 32768 : psCtor_p->correctX1D(correction, iy+inner/2);
149 32768 : Int k=0;
150 16809984 : for (Int ix=-inner/2;ix<inner/2;ix++)
151 : {
152 16777216 : screen(ix+convSize/2,iy+convSize/2)=correction(k++);//correction(ix+inner/2);
153 : }
154 : }
155 : }
156 : (void)sampling;
157 64 : }
158 : //
159 : //-------------------------------------------------------------------------
160 : //
161 0 : void PSTerm::normalizeImage(Lattice<Complex>& skyImage,
162 : const Matrix<Float>& weights)
163 : {
164 0 : Int inx = skyImage.shape()(0);
165 : // Int iny = skyImage.shape()(1);
166 0 : Vector<Complex> correction(inx);
167 0 : correction=Complex(1.0, 0.0);
168 0 : IPosition cursorShape(4, inx, 1, 1, 1);
169 0 : IPosition axisPath(4, 0, 1, 2, 3);
170 0 : LatticeStepper lsx(skyImage.shape(), cursorShape, axisPath);
171 0 : LatticeIterator<Complex> lix(skyImage, lsx);
172 :
173 0 : for(lix.reset();!lix.atEnd();lix++)
174 : {
175 0 : Int pol=lix.position()(2), chan=lix.position()(3);
176 :
177 0 : if(weights(pol, chan)!=0.0)
178 : {
179 0 : Int iy=lix.position()(1);
180 0 : psCtor_p->correctX1D(correction,iy);
181 : // for (Int ix=0;ix<inx;ix++)
182 : // correction(ix)*=sincConv(ix)*sincConv(iy);
183 :
184 0 : lix.rwVectorCursor()/=correction;
185 : }
186 : else
187 0 : lix.woCursor()=0.0;
188 : }
189 0 : }
190 : };
|
