Line data Source code
1 : //# SDDataSampling.cc: Implementation of SDDataSampling class
2 : //# Copyright (C) 1997,1998,1999,2000,2001,2003
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 addressed 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 : //# $Id$
27 :
28 : #include <synthesis/DataSampling/SDDataSampling.h>
29 : #include <casacore/coordinates/Coordinates/CoordinateSystem.h>
30 : #include <casacore/coordinates/Coordinates/DirectionCoordinate.h>
31 : #include <casacore/ms/MeasurementSets/MSColumns.h>
32 : #include <casacore/casa/BasicSL/Constants.h>
33 : #include <synthesis/TransformMachines/SkyJones.h>
34 : #include <msvis/MSVis/VisSet.h>
35 : #include <msvis/MSVis/VisBuffer.h>
36 : #include <msvis/MSVis/VisibilityIterator.h>
37 : #include <casacore/images/Images/ImageInterface.h>
38 : #include <casacore/images/Images/PagedImage.h>
39 : #include <casacore/images/Images/TempImage.h>
40 : #include <casacore/casa/Arrays/ArrayLogical.h>
41 : #include <casacore/casa/Arrays/ArrayMath.h>
42 : #include <casacore/casa/Arrays/MaskedArray.h>
43 : #include <casacore/casa/Arrays/Array.h>
44 : #include <casacore/casa/Arrays/Array.h>
45 : #include <casacore/casa/Arrays/Vector.h>
46 : #include <casacore/casa/Arrays/Matrix.h>
47 : #include <casacore/casa/BasicSL/String.h>
48 : #include <casacore/casa/Utilities/Assert.h>
49 : #include <casacore/casa/Exceptions/Error.h>
50 : #include <casacore/casa/System/ProgressMeter.h>
51 : #include <sstream>
52 :
53 : using namespace casacore;
54 : namespace casa { //# NAMESPACE CASA - BEGIN
55 :
56 0 : SDDataSampling::SDDataSampling(MeasurementSet& ms,
57 : SkyJones& sj,
58 : const CoordinateSystem& coords,
59 : const IPosition& shape,
60 0 : const Quantity& sigmaConst)
61 : {
62 :
63 0 : LogIO os(LogOrigin("SDDataSampling", "SDDataSampling"));
64 :
65 0 : DataSampling::IDLScript_p="@app_sd";
66 :
67 : // Now create the VisSet
68 0 : Block<int> sort(1);
69 0 : sort[0] = MS::TIME;
70 :
71 0 : Matrix<Int> noselection;
72 0 : VisSet vs(ms,sort,noselection);
73 :
74 : // First get the CoordinateSystem for the image and then find
75 : // the DirectionCoordinate
76 0 : Int directionIndex=coords.findCoordinate(Coordinate::DIRECTION);
77 0 : AlwaysAssert(directionIndex>=0, AipsError);
78 0 : DirectionCoordinate directionCoord=coords.directionCoordinate(directionIndex);
79 :
80 0 : Int lastRow = 0;
81 :
82 0 : MSPointingColumns mspc(ms.pointing());
83 :
84 0 : lastIndex_p=0;
85 :
86 0 : VisIter& vi=vs.iter();
87 :
88 : // Get bigger chunks o'data: this should be tuned some time
89 : // since it may be wrong for e.g. spectral line
90 0 : vi.setRowBlocking(100);
91 :
92 0 : VisBuffer vb(vi);
93 0 : vi.originChunks();
94 0 : vi.origin();
95 :
96 0 : Vector<Double> imagePos;
97 0 : imagePos = directionCoord.referenceValue();
98 :
99 0 : Vector<Double> worldPos;
100 :
101 0 : MDirection imagePosMeas;
102 0 : MDirection worldPosMeas;
103 :
104 : // First try the POINTING sub-table
105 0 : Int pointIndex=getIndex(mspc, vb.time()(0));
106 : // If no valid POINTING entry, then use FIELD phase center
107 0 : MSColumns msc(ms);
108 0 : if(pointIndex >= 0 || pointIndex < static_cast<Int>(mspc.time().nrow()))
109 0 : worldPosMeas = mspc.directionMeas(pointIndex);
110 : else
111 0 : worldPosMeas = msc.field().phaseDirMeas(vb.fieldId());
112 :
113 : MDirection::Convert pointingToImage(worldPosMeas,
114 0 : directionCoord.directionType());
115 0 : imagePosMeas = pointingToImage(worldPosMeas);
116 0 : directionCoord.setReferenceValue(imagePosMeas.getAngle().getValue());
117 :
118 0 : Vector<Double> unitVec(2);
119 0 : Int nx=shape(0);
120 0 : Int ny=shape(1);
121 0 : unitVec(0)=nx/2;
122 0 : unitVec(1)=ny/2;
123 0 : directionCoord.setReferencePixel(unitVec);
124 0 : CoordinateSystem iCoords(coords);
125 0 : iCoords.replaceCoordinate(directionCoord, directionIndex);
126 :
127 0 : TempImage<Float> PB(shape, iCoords);
128 0 : PB.set(1.0);
129 :
130 0 : sj.applySquare(PB, PB, vb, 0);
131 0 : prf_p=PB.getSlice(IPosition(4, 0, 0, 0, 0),
132 0 : IPosition(4, nx, ny, 1, 1), true);
133 :
134 : // Reset the direction coordinate
135 0 : directionCoord=coords.directionCoordinate(directionIndex);
136 : // Now fill in the data and position columns
137 0 : ProgressMeter pm(1.0, Double(ms.nrow()), "Sampling Data", "", "", "", true);
138 :
139 : // Loop over all visibilities
140 0 : Int cohDone = 0;
141 0 : Float sigmaVal=sigmaConst.getValue();
142 0 : Vector<Double> xyPos(2);
143 :
144 0 : Array<Float> dx(IPosition(2, 2, ms.nrow()));
145 0 : dx=-1.0;
146 0 : Array<Float> data(IPosition(1, ms.nrow()));
147 0 : data=0.0;
148 0 : Array<Float> sigma(IPosition(1, ms.nrow()));
149 0 : sigma=-1.0;
150 :
151 0 : for (vi.originChunks();vi.moreChunks();vi.nextChunk()) {
152 0 : for (vi.origin(); vi.more(); vi++) {
153 0 : for (Int row=0;row<vb.nRow();row++) {
154 0 : Int pointIndex=getIndex(mspc, vb.time()(row));
155 0 : if(pointIndex >= 0 || pointIndex < static_cast<Int>(mspc.time().nrow())){
156 : imagePosMeas =
157 0 : pointingToImage(mspc.directionMeas(pointIndex));
158 0 : if(directionCoord.toPixel(xyPos, imagePosMeas)) {
159 0 : if((!vb.flagRow()(row))&&
160 0 : (vb.sigma()(row)>0.0)&&
161 0 : (Float(xyPos(0))>0.0)&&
162 0 : (Float(xyPos(1))>0.0)) {
163 0 : dx(IPosition(2, 0, lastRow)) = Float(xyPos(0));
164 0 : dx(IPosition(2, 1, lastRow)) = Float(xyPos(1));
165 0 : IPosition irow(1, lastRow);
166 0 : data(irow) = real(vb.correctedVisCube()(0,0,row));
167 0 : if(sigmaVal>0.0) {
168 0 : sigma(irow) = sigmaVal;
169 : }
170 : else {
171 0 : sigma(irow) = vb.sigma()(row);
172 : }
173 0 : lastRow++;
174 : }
175 : }
176 : }
177 : }
178 0 : cohDone+=vb.nRow();
179 0 : pm.update(Double(cohDone));
180 : }
181 : }
182 0 : if(lastRow==0) {
183 0 : LogIO os(LogOrigin("SDDataSampling", "SDDataSampling()", WHERE));
184 : os << LogIO::SEVERE << "No valid data: check image parameters, sigmas in data"
185 0 : << LogIO::EXCEPTION;
186 : }
187 : // Now copy good rows to output arrays
188 0 : dx_p.resize(IPosition(2, 2, lastRow));
189 0 : data_p.resize(IPosition(1, lastRow));
190 0 : sigma_p.resize(IPosition(1, lastRow));
191 0 : for (Int row=0;row<lastRow;row++) {
192 0 : IPosition ip(2, 0, row);
193 0 : dx_p(ip)=dx(ip);
194 0 : ip(0)=1;
195 0 : dx_p(ip)=dx(ip);
196 0 : IPosition rp(1, row);
197 0 : sigma_p(rp)=sigma(rp);
198 0 : data_p(rp)=data(rp);
199 : }
200 :
201 0 : }
202 :
203 : //----------------------------------------------------------------------
204 0 : SDDataSampling& SDDataSampling::operator=(const SDDataSampling& other)
205 : {
206 : if(this!=&other) {
207 : };
208 0 : return *this;
209 : };
210 :
211 : //----------------------------------------------------------------------
212 0 : SDDataSampling::SDDataSampling(const SDDataSampling& other)
213 0 : :DataSampling(other)
214 : {
215 0 : operator=(other);
216 0 : }
217 :
218 : //----------------------------------------------------------------------
219 0 : SDDataSampling::~SDDataSampling() {
220 0 : }
221 :
222 0 : void SDDataSampling::ok() {
223 0 : }
224 :
225 0 : Int SDDataSampling::getIndex(const MSPointingColumns& mspc, const Double& time) {
226 0 : Int start=lastIndex_p;
227 : // Search forwards
228 0 : Int nrows=mspc.time().nrow();
229 0 : for (Int i=start;i<nrows;i++) {
230 0 : if(abs(mspc.time()(i)-time) < mspc.interval()(i)) {
231 0 : lastIndex_p=i;
232 0 : return i;
233 : }
234 : }
235 : // Search backwards
236 0 : for (Int i=start;i>=0;i--) {
237 0 : if(abs(mspc.time()(i)-time) < mspc.interval()(i)) {
238 0 : lastIndex_p=i;
239 0 : return i;
240 : }
241 : }
242 : // No match!
243 0 : return -1;
244 : }
245 :
246 : } //# NAMESPACE CASA - END
247 :
|
