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