Commits
65 65 | |
66 66 | |
67 67 | |
68 68 | |
69 69 | |
70 70 | |
71 71 | |
72 72 | |
73 73 | |
74 74 | using namespace casacore; |
75 - | namespace casa { //# NAMESPACE CASA - BEGIN |
76 - | |
77 - | // Public functions |
75 + | namespace casa { |
78 76 | |
79 77 | ImagePolarimetry::ImagePolarimetry (const ImageInterface<Float>& image) |
80 - | : itsFitterPtr(0), |
81 - | itsOldClip(0.0), |
82 - | _beamsEqMat() |
78 + | : _image(image.cloneII()) |
83 79 | { |
84 - | itsInImagePtr.reset(image.cloneII()); |
85 - | // |
86 - | itsStokesPtr.resize(4); |
87 - | itsStokesStatsPtr.resize(4); |
88 - | itsStokesPtr.set(0); |
89 - | itsStokesStatsPtr.set(0); |
90 - | findStokes(); |
80 + | _stokes.resize(4); |
81 + | _stokesStats.resize(4); |
82 + | _stokes.set(0); |
83 + | _stokesStats.set(0); |
84 + | _findStokes(); |
91 85 | _createBeamsEqMat(); |
92 86 | } |
93 87 | |
94 88 | |
95 89 | ImagePolarimetry::ImagePolarimetry(const ImagePolarimetry &other) |
96 90 | { |
97 91 | operator=(other); |
98 92 | } |
99 93 | |
100 94 | |
101 95 | ImagePolarimetry &ImagePolarimetry::operator=(const ImagePolarimetry &other) |
102 96 | { |
103 97 | if (this != &other) { |
104 - | itsInImagePtr.reset(other.itsInImagePtr->cloneII()); |
98 + | _image.reset(other._image->cloneII()); |
105 99 | // |
106 - | const uInt n = itsStokesPtr.nelements(); |
100 + | const uInt n = _stokes.nelements(); |
107 101 | for (uInt i=0; i<n; i++) { |
108 - | if (itsStokesPtr[i]!=0) { |
109 - | delete itsStokesPtr[i]; |
110 - | itsStokesPtr[i] = 0; |
102 + | if (_stokes[i]!=0) { |
103 + | delete _stokes[i]; |
104 + | _stokes[i] = 0; |
111 105 | } |
112 - | if (other.itsStokesPtr[i]!=0) { |
113 - | itsStokesPtr[i] = other.itsStokesPtr[i]->cloneII(); |
106 + | if (other._stokes[i]!=0) { |
107 + | _stokes[i] = other._stokes[i]->cloneII(); |
114 108 | } |
115 109 | } |
116 110 | |
117 111 | // Just delete fitter. It will make a new one when needed. |
118 112 | |
119 - | if (itsFitterPtr!= 0) { |
120 - | delete itsFitterPtr; |
121 - | itsFitterPtr = 0; |
113 + | if (_fitter!= 0) { |
114 + | delete _fitter; |
115 + | _fitter = 0; |
122 116 | } |
123 117 | |
124 118 | // Remake Statistics objects as needed |
125 119 | |
126 - | itsOldClip = 0.0; |
120 + | _oldClip = 0.0; |
127 121 | for (uInt i=0; i<n; i++) { |
128 - | if (itsStokesStatsPtr[i]!=0) { |
129 - | delete itsStokesStatsPtr[i]; |
130 - | itsStokesStatsPtr[i] = 0; |
122 + | if (_stokesStats[i]!=0) { |
123 + | delete _stokesStats[i]; |
124 + | _stokesStats[i] = 0; |
131 125 | } |
132 126 | } |
133 127 | _beamsEqMat.assign(other._beamsEqMat); |
134 128 | } |
135 129 | return *this; |
136 130 | } |
137 131 | |
138 132 | |
139 133 | ImagePolarimetry::~ImagePolarimetry() |
140 134 | { |
141 - | cleanup(); |
135 + | _cleanup(); |
142 136 | } |
143 137 | |
144 138 | |
145 139 | // Public methods |
146 140 | |
147 141 | ImageExpr<Complex> ImagePolarimetry::complexLinearPolarization() |
148 142 | { |
149 - | hasQU(); |
143 + | _hasQU(); |
150 144 | _checkQUBeams(false); |
151 - | LatticeExprNode node(formComplex(*itsStokesPtr[ImagePolarimetry::Q], |
152 - | *itsStokesPtr[ImagePolarimetry::U])); |
145 + | LatticeExprNode node(formComplex(*_stokes[ImagePolarimetry::Q], |
146 + | *_stokes[ImagePolarimetry::U])); |
153 147 | LatticeExpr<Complex> le(node); |
154 148 | ImageExpr<Complex> ie(le, String("ComplexLinearPolarization")); |
155 - | fiddleStokesCoordinate(ie, Stokes::Plinear); // Need a Complex Linear Polarization type |
149 + | _fiddleStokesCoordinate(ie, Stokes::Plinear); // Need a Complex Linear Polarization type |
156 150 | // |
157 - | ie.setUnits(itsInImagePtr->units()); |
151 + | ie.setUnits(_image->units()); |
158 152 | _setInfo(ie, Q); |
159 153 | return ie; |
160 154 | } |
161 155 | |
162 156 | void ImagePolarimetry::_setInfo(ImageInterface<Complex>& im, const StokesTypes stokes) const { |
163 - | ImageInfo info = itsInImagePtr->imageInfo(); |
157 + | ImageInfo info = _image->imageInfo(); |
164 158 | if (info.hasMultipleBeams()) { |
165 - | info.setBeams(itsStokesPtr[stokes]->imageInfo().getBeamSet()); |
159 + | info.setBeams(_stokes[stokes]->imageInfo().getBeamSet()); |
166 160 | } |
167 161 | im.setImageInfo(info); |
168 162 | } |
169 163 | |
170 164 | void ImagePolarimetry::_setInfo(ImageInterface<Float>& im, const StokesTypes stokes) const { |
171 - | ImageInfo info = itsInImagePtr->imageInfo(); |
165 + | ImageInfo info = _image->imageInfo(); |
172 166 | if (info.hasMultipleBeams()) { |
173 - | info.setBeams(itsStokesPtr[stokes]->imageInfo().getBeamSet()); |
167 + | info.setBeams(_stokes[stokes]->imageInfo().getBeamSet()); |
174 168 | } |
175 169 | im.setImageInfo(info); |
176 170 | } |
177 171 | |
178 172 | |
179 173 | ImageExpr<Complex> ImagePolarimetry::complexFractionalLinearPolarization() |
180 174 | { |
181 175 | LogIO os(LogOrigin("ImagePolarimetry", __func__, WHERE)); |
182 - | hasQU(); |
183 - | if (itsStokesPtr[ImagePolarimetry::I]==0) { |
176 + | _hasQU(); |
177 + | if (_stokes[ImagePolarimetry::I]==0) { |
184 178 | os << "This image does not have Stokes I so cannot provide fractional linear polarization" << LogIO::EXCEPTION; |
185 179 | } |
186 180 | |
187 181 | _checkIQUBeams(false); |
188 - | LatticeExprNode nodeQU(formComplex(*itsStokesPtr[ImagePolarimetry::Q], |
189 - | *itsStokesPtr[ImagePolarimetry::U])); |
190 - | LatticeExprNode nodeI(*itsStokesPtr[ImagePolarimetry::I]); |
182 + | LatticeExprNode nodeQU(formComplex(*_stokes[ImagePolarimetry::Q], |
183 + | *_stokes[ImagePolarimetry::U])); |
184 + | LatticeExprNode nodeI(*_stokes[ImagePolarimetry::I]); |
191 185 | LatticeExpr<Complex> le(nodeQU/nodeI); |
192 186 | ImageExpr<Complex> ie(le, String("ComplexFractionalLinearPolarization")); |
193 - | fiddleStokesCoordinate(ie, Stokes::PFlinear); // Need a Complex Linear Polarization type |
187 + | _fiddleStokesCoordinate(ie, Stokes::PFlinear); // Need a Complex Linear Polarization type |
194 188 | ie.setUnits(Unit("")); |
195 189 | _setInfo(ie, I); |
196 190 | return ie; |
197 191 | } |
198 192 | |
199 193 | ImageExpr<Float> ImagePolarimetry::fracLinPol(Bool debias, Float clip, Float sigma) |
200 194 | { |
201 195 | LogIO os(LogOrigin("ImagePolarimetry", __func__, WHERE)); |
202 - | hasQU(); |
196 + | _hasQU(); |
203 197 | |
204 - | if (itsStokesPtr[ImagePolarimetry::I]==0) { |
198 + | if (_stokes[ImagePolarimetry::I]==0) { |
205 199 | os << "This image does not have Stokes I so cannot provide fractional linear polarization" << LogIO::EXCEPTION; |
206 200 | } |
207 201 | Vector<StokesTypes> types(3); |
208 202 | types[0] = I; types[1] = Q; types[2] = U; |
209 203 | _checkIQUBeams(false); |
210 204 | // Make nodes |
211 205 | |
212 - | LatticeExprNode nodePol = makePolIntNode(os, debias, clip, sigma, true, false); |
213 - | LatticeExprNode nodeI(*itsStokesPtr[ImagePolarimetry::I]); |
206 + | LatticeExprNode nodePol = _makePolIntNode(os, debias, clip, sigma, true, false); |
207 + | LatticeExprNode nodeI(*_stokes[ImagePolarimetry::I]); |
214 208 | |
215 209 | // Make expression |
216 210 | |
217 211 | LatticeExpr<Float> le(nodePol/nodeI); |
218 212 | ImageExpr<Float> ie(le, String("FractionalLinearPolarization")); |
219 213 | // |
220 214 | ie.setUnits(Unit("")); |
221 - | ImageInfo ii = itsInImagePtr->imageInfo(); |
215 + | ImageInfo ii = _image->imageInfo(); |
222 216 | ii.removeRestoringBeam(); |
223 217 | ie.setImageInfo(ii); |
224 - | |
225 - | // Fiddle Stokes coordinate in ImageExpr |
226 - | |
227 - | fiddleStokesCoordinate(ie, Stokes::PFlinear); |
218 + | _fiddleStokesCoordinate(ie, Stokes::PFlinear); |
228 219 | // |
229 220 | return ie; |
230 221 | } |
231 222 | |
232 223 | |
233 224 | ImageExpr<Float> ImagePolarimetry::sigmaFracLinPol(Float clip, Float sigma) |
234 225 | // |
235 226 | // sigma_m = m * sqrt( (sigmaP/p)**2 + (sigmaI/I)**2) ) |
236 227 | // sigmaP = sigmaQU |
237 228 | // sigmaI = sigmaI |
238 229 | // |
239 230 | { |
240 231 | LogIO os(LogOrigin("ImagePolarimetry", __func__, WHERE)); |
241 - | hasQU(); |
242 - | if (itsStokesPtr[ImagePolarimetry::I]==0) { |
232 + | _hasQU(); |
233 + | if (_stokes[ImagePolarimetry::I]==0) { |
243 234 | os << "This image does not have Stokes I so cannot provide fractional linear polarization" << LogIO::EXCEPTION; |
244 235 | } |
245 236 | _checkIQUBeams(false); |
246 237 | // Make nodes. Don't bother debiasing. |
247 238 | |
248 239 | Bool debias = false; |
249 - | LatticeExprNode nodePol = makePolIntNode(os, debias, clip, sigma, true, false); |
250 - | LatticeExprNode nodeI(*itsStokesPtr[ImagePolarimetry::I]); |
240 + | LatticeExprNode nodePol = _makePolIntNode(os, debias, clip, sigma, true, false); |
241 + | LatticeExprNode nodeI(*_stokes[ImagePolarimetry::I]); |
251 242 | |
252 243 | // Make expression. We assume sigmaI = sigmaQU which is true with |
253 244 | // no dynamic range limititation. Perhaps we should work out |
254 245 | // sigmaI as well. |
255 246 | |
256 247 | Float sigma2 = sigmaLinPolInt(clip, sigma); |
257 248 | LatticeExprNode n0(nodePol / nodeI); |
258 249 | LatticeExprNode n1(pow(sigma2/nodePol,2)); |
259 250 | LatticeExprNode n2(pow(sigma2/nodeI,2)); |
260 251 | LatticeExpr<Float> le(n0 * sqrt(n1 + n2)); |
261 252 | ImageExpr<Float> ie(le, String("FractionalLinearPolarizationError")); |
262 253 | // |
263 254 | ie.setUnits(Unit("")); |
264 - | ImageInfo ii = itsInImagePtr->imageInfo(); |
255 + | ImageInfo ii = _image->imageInfo(); |
265 256 | ii.removeRestoringBeam(); |
266 257 | ie.setImageInfo(ii); |
267 - | |
268 - | // Fiddle Stokes coordinate in ImageExpr |
269 - | |
270 - | fiddleStokesCoordinate(ie, Stokes::PFlinear); |
271 - | // |
258 + | _fiddleStokesCoordinate(ie, Stokes::PFlinear); |
272 259 | return ie; |
273 260 | } |
274 261 | |
275 262 | |
276 263 | ImageExpr<Float> ImagePolarimetry::fracTotPol(Bool debias, Float clip, Float sigma) |
277 264 | { |
278 265 | LogIO os(LogOrigin("ImagePolarimetry", __func__, WHERE)); |
279 266 | Bool doLin, doCirc; |
280 267 | _setDoLinDoCirc(doLin, doCirc); |
281 268 | |
282 269 | // Make nodes |
283 270 | |
284 - | LatticeExprNode nodePol = makePolIntNode(os, debias, clip, sigma, doLin, doCirc); |
285 - | LatticeExprNode nodeI(*itsStokesPtr[ImagePolarimetry::I]); |
271 + | LatticeExprNode nodePol = _makePolIntNode(os, debias, clip, sigma, doLin, doCirc); |
272 + | LatticeExprNode nodeI(*_stokes[ImagePolarimetry::I]); |
286 273 | |
287 274 | // Make expression |
288 275 | |
289 276 | LatticeExpr<Float> le(nodePol/nodeI); |
290 277 | ImageExpr<Float> ie(le, String("FractionalTotalPolarization")); |
291 278 | // |
292 279 | ie.setUnits(Unit("")); |
293 - | ImageInfo ii = itsInImagePtr->imageInfo(); |
280 + | ImageInfo ii = _image->imageInfo(); |
294 281 | ii.removeRestoringBeam(); |
295 282 | ie.setImageInfo(ii); |
296 - | |
297 - | // Fiddle Stokes coordinate in ImageExpr |
298 - | |
299 - | fiddleStokesCoordinate(ie, Stokes::PFtotal); |
300 - | // |
283 + | _fiddleStokesCoordinate(ie, Stokes::PFtotal); |
301 284 | return ie; |
302 285 | } |
303 286 | |
304 287 | void ImagePolarimetry::_setDoLinDoCirc(Bool& doLin, Bool& doCirc) const { |
305 288 | LogIO os(LogOrigin("ImagePolarimetry", __func__, WHERE)); |
306 289 | doLin = ( |
307 - | itsStokesPtr[ImagePolarimetry::Q]!=0 |
308 - | && itsStokesPtr[ImagePolarimetry::U]!=0 |
290 + | _stokes[ImagePolarimetry::Q]!=0 |
291 + | && _stokes[ImagePolarimetry::U]!=0 |
309 292 | ); |
310 - | doCirc = (itsStokesPtr[ImagePolarimetry::V]!=0); |
293 + | doCirc = (_stokes[ImagePolarimetry::V]!=0); |
311 294 | AlwaysAssert((doLin||doCirc), AipsError); // Should never happen |
312 - | if (itsStokesPtr[ImagePolarimetry::I]==0) { |
295 + | if (_stokes[ImagePolarimetry::I]==0) { |
313 296 | os << "This image does not have Stokes I so this calculation cannot be carried out" |
314 297 | << LogIO::EXCEPTION; |
315 298 | } |
316 299 | if (doLin) { |
317 300 | if (! _checkIQUBeams(false, false)) { |
318 301 | os << LogIO::WARN |
319 302 | << "I, Q, and U beams are not the same, cannot do linear portion" |
320 303 | << LogIO::POST; |
321 304 | doLin = false; |
322 305 | } |
342 325 | // sigmaI = sigmaI |
343 326 | // |
344 327 | { |
345 328 | LogIO os(LogOrigin("ImagePolarimetry", __func__, WHERE)); |
346 329 | Bool doLin, doCirc; |
347 330 | _setDoLinDoCirc(doLin, doCirc); |
348 331 | |
349 332 | // Make nodes. Don't bother debiasing. |
350 333 | |
351 334 | Bool debias = false; |
352 - | LatticeExprNode nodePol = makePolIntNode(os, debias, clip, sigma, doLin, doCirc); |
353 - | LatticeExprNode nodeI(*itsStokesPtr[ImagePolarimetry::I]); |
335 + | LatticeExprNode nodePol = _makePolIntNode(os, debias, clip, sigma, doLin, doCirc); |
336 + | LatticeExprNode nodeI(*_stokes[ImagePolarimetry::I]); |
354 337 | |
355 338 | // Make expression. We assume sigmaI = sigmaQU which is true with |
356 339 | // no dynamic range limitation. Perhaps we should work out |
357 340 | // sigmaI as well. |
358 341 | |
359 342 | Float sigma2 = sigmaTotPolInt(clip, sigma); |
360 343 | LatticeExprNode n0(nodePol / nodeI); |
361 344 | LatticeExprNode n1(pow(sigma2/nodePol,2)); |
362 345 | LatticeExprNode n2(pow(sigma2/nodeI,2)); |
363 346 | LatticeExpr<Float> le(n0 * sqrt(n1 + n2)); |
364 347 | ImageExpr<Float> ie(le, String("FractionalLinearPolarizationError")); |
365 348 | // |
366 349 | ie.setUnits(Unit("")); |
367 - | ImageInfo ii = itsInImagePtr->imageInfo(); |
350 + | ImageInfo ii = _image->imageInfo(); |
368 351 | ii.removeRestoringBeam(); |
369 352 | ie.setImageInfo(ii); |
370 - | |
371 - | // Fiddle Stokes coordinate in ImageExpr |
372 - | |
373 - | fiddleStokesCoordinate(ie, Stokes::PFlinear); |
374 - | // |
353 + | _fiddleStokesCoordinate(ie, Stokes::PFlinear); |
375 354 | return ie; |
376 355 | } |
377 356 | |
378 357 | |
379 358 | |
380 359 | void ImagePolarimetry::fourierRotationMeasure(ImageInterface<Complex>& cpol, |
381 360 | Bool zeroZeroLag) |
382 361 | { |
383 362 | LogIO os(LogOrigin("ImagePolarimetry", __func__, WHERE)); |
384 - | hasQU(); |
363 + | _hasQU(); |
385 364 | _checkQUBeams(true, false); |
386 365 | // Check image shape |
387 366 | CoordinateSystem dCS; |
388 367 | Stokes::StokesTypes dType = Stokes::Plinear; |
389 368 | IPosition shape = singleStokesShape(dCS, dType); |
390 369 | if (!cpol.shape().isEqual(shape)) { |
391 370 | os << "The provided image has the wrong shape " << cpol.shape() << endl; |
392 371 | os << "It should be of shape " << shape << LogIO::EXCEPTION; |
393 372 | } |
394 373 | // Find spectral coordinate. |
395 374 | |
396 - | const CoordinateSystem& cSys = itsInImagePtr->coordinates(); |
375 + | const CoordinateSystem& cSys = _image->coordinates(); |
397 376 | Int spectralCoord, fAxis; |
398 - | findFrequencyAxis (spectralCoord, fAxis, cSys, -1); |
377 + | _findFrequencyAxis (spectralCoord, fAxis, cSys, -1); |
399 378 | |
400 379 | // Make Complex (Q,U) image |
401 380 | LatticeExprNode node; |
402 381 | if (zeroZeroLag) { |
403 - | TempImage<Float> tQ(itsStokesPtr[ImagePolarimetry::Q]->shape(), |
404 - | itsStokesPtr[ImagePolarimetry::Q]->coordinates()); |
405 - | if (itsStokesPtr[ImagePolarimetry::Q]->isMasked()) { |
382 + | TempImage<Float> tQ(_stokes[ImagePolarimetry::Q]->shape(), |
383 + | _stokes[ImagePolarimetry::Q]->coordinates()); |
384 + | if (_stokes[ImagePolarimetry::Q]->isMasked()) { |
406 385 | tQ.makeMask(String("mask0"), true, true, false, false); |
407 386 | } |
408 - | copyDataAndMask (tQ, *itsStokesPtr[ImagePolarimetry::Q]); |
409 - | subtractProfileMean (tQ, fAxis); |
387 + | _copyDataAndMask (tQ, *_stokes[ImagePolarimetry::Q]); |
388 + | _subtractProfileMean (tQ, fAxis); |
410 389 | // |
411 - | TempImage<Float> tU(itsStokesPtr[ImagePolarimetry::U]->shape(), |
412 - | itsStokesPtr[ImagePolarimetry::U]->coordinates()); |
413 - | if (itsStokesPtr[ImagePolarimetry::U]->isMasked()) { |
390 + | TempImage<Float> tU(_stokes[ImagePolarimetry::U]->shape(), |
391 + | _stokes[ImagePolarimetry::U]->coordinates()); |
392 + | if (_stokes[ImagePolarimetry::U]->isMasked()) { |
414 393 | tU.makeMask(String("mask0"), true, true, false, false); |
415 394 | } |
416 - | copyDataAndMask (tU, *itsStokesPtr[ImagePolarimetry::U]); |
417 - | subtractProfileMean (tU, fAxis); |
395 + | _copyDataAndMask (tU, *_stokes[ImagePolarimetry::U]); |
396 + | _subtractProfileMean (tU, fAxis); |
418 397 | |
419 398 | // The TempImages will be cloned be LatticeExprNode so it's ok |
420 399 | // that they go out of scope |
421 400 | |
422 401 | node = LatticeExprNode(formComplex(tQ, tU)); |
423 402 | } else { |
424 - | node = LatticeExprNode(formComplex(*itsStokesPtr[ImagePolarimetry::Q], |
425 - | *itsStokesPtr[ImagePolarimetry::U])); |
403 + | node = LatticeExprNode(formComplex(*_stokes[ImagePolarimetry::Q], |
404 + | *_stokes[ImagePolarimetry::U])); |
426 405 | } |
427 406 | LatticeExpr<Complex> le(node); |
428 407 | ImageExpr<Complex> ie(le, String("ComplexLinearPolarization")); |
429 408 | // Do FFT of spectral coordinate |
430 409 | |
431 410 | Vector<Bool> axes(ie.ndim(),false); |
432 411 | axes(fAxis) = true; |
433 412 | ImageFFT<Complex> fftserver; |
434 413 | fftserver.fft(ie, axes); |
435 414 | // Recover Complex result. Coordinates are updated to include Fourier coordinate, |
436 415 | // miscellaneous things (MiscInfo, ImageInfo, units, history) and mask |
437 416 | // (if output has one) are copied to cpol |
438 417 | cout << "*** names " << cpol.coordinates().worldAxisNames() << endl; |
439 418 | |
440 419 | fftserver.getComplex(cpol); |
441 420 | // Fiddle time coordinate to be a RotationMeasure coordinate |
442 421 | |
443 - | Quantum<Double> f = findCentralFrequency(cSys.coordinate(spectralCoord), ie.shape()(fAxis)); |
444 - | fiddleTimeCoordinate(cpol, f, spectralCoord); |
445 - | |
446 - | // Set Stokes coordinate to be correct type |
447 - | fiddleStokesCoordinate(cpol, Stokes::Plinear); |
422 + | Quantum<Double> f = _findCentralFrequency(cSys.coordinate(spectralCoord), ie.shape()(fAxis)); |
423 + | _fiddleTimeCoordinate(cpol, f, spectralCoord); |
424 + | // Set Stokes coordinate to be correct type |
425 + | _fiddleStokesCoordinate(cpol, Stokes::Plinear); |
448 426 | |
449 427 | // Set units and ImageInfo |
450 428 | |
451 - | cpol.setUnits(itsInImagePtr->units()); |
429 + | cpol.setUnits(_image->units()); |
452 430 | _setInfo(cpol, Q); |
453 431 | } |
454 432 | |
455 433 | Float ImagePolarimetry::sigmaLinPolInt(Float clip, Float sigma) |
456 434 | // |
457 435 | // sigma_P = sigma_QU |
458 436 | // |
459 437 | { |
460 438 | LogIO os(LogOrigin("ImagePolarimetry", __func__, WHERE)); |
461 - | if (itsStokesPtr[ImagePolarimetry::Q]==0 && itsStokesPtr[ImagePolarimetry::U]==0) { |
439 + | if (_stokes[ImagePolarimetry::Q]==0 && _stokes[ImagePolarimetry::U]==0) { |
462 440 | os << "This image does not have Stokes Q and U so cannot provide linear polarization" << LogIO::EXCEPTION; |
463 441 | } |
464 442 | _checkQUBeams(false); |
465 443 | // Get value |
466 444 | |
467 445 | Float sigma2 = 0.0; |
468 446 | if (sigma > 0) { |
469 447 | sigma2 = sigma; |
470 448 | } else { |
471 449 | os << LogIO::NORMAL << "Determined noise from Q&U images to be "; |
472 - | Float sq = ImagePolarimetry::sigma(ImagePolarimetry::Q, clip); |
473 - | Float su = ImagePolarimetry::sigma(ImagePolarimetry::U, clip); |
450 + | Float sq = _sigma(ImagePolarimetry::Q, clip); |
451 + | Float su = _sigma(ImagePolarimetry::U, clip); |
474 452 | sigma2 = (sq+su)/2.0; |
475 453 | } |
476 454 | return sigma2; |
477 455 | } |
478 456 | |
479 457 | ImageExpr<Float> ImagePolarimetry::linPolPosAng(Bool radians) const |
480 458 | { |
481 459 | LogIO os(LogOrigin("ImagePolarimetry", __func__, WHERE)); |
482 - | if (itsStokesPtr[ImagePolarimetry::Q]==0 && itsStokesPtr[ImagePolarimetry::U]==0) { |
460 + | if (_stokes[ImagePolarimetry::Q]==0 && _stokes[ImagePolarimetry::U]==0) { |
483 461 | os << "This image does not have Stokes Q and U so cannot provide linear polarization" << LogIO::EXCEPTION; |
484 462 | } |
485 463 | _checkQUBeams(false); |
486 464 | |
487 465 | // Make expression. LEL function "pa" returns degrees |
488 466 | |
489 467 | Float fac = 1.0; |
490 468 | if (radians) fac = C::pi / 180.0; |
491 - | LatticeExprNode node(fac*pa(*itsStokesPtr[ImagePolarimetry::U], |
492 - | *itsStokesPtr[ImagePolarimetry::Q])); |
469 + | LatticeExprNode node(fac*pa(*_stokes[ImagePolarimetry::U], |
470 + | *_stokes[ImagePolarimetry::Q])); |
493 471 | LatticeExpr<Float> le(node); |
494 472 | ImageExpr<Float> ie(le, String("LinearlyPolarizedPositionAngle")); |
495 473 | // |
496 474 | if (radians) { |
497 475 | ie.setUnits(Unit("rad")); |
498 476 | } else { |
499 477 | ie.setUnits(Unit("deg")); |
500 478 | } |
501 - | ImageInfo ii = itsInImagePtr->imageInfo(); |
479 + | ImageInfo ii = _image->imageInfo(); |
502 480 | ii.removeRestoringBeam(); |
503 481 | ie.setImageInfo(ii); |
504 - | |
505 - | // Fiddle Stokes coordinate |
506 - | |
507 - | fiddleStokesCoordinate(ie, Stokes::Pangle); |
508 - | |
482 + | _fiddleStokesCoordinate(ie, Stokes::Pangle); |
509 483 | return ie; |
510 484 | } |
511 485 | |
512 486 | ImageExpr<Float> ImagePolarimetry::sigmaLinPolPosAng(Bool radians, Float clip, Float sigma) |
513 487 | // |
514 488 | // sigma_PA = sigmaQU / 2P |
515 489 | // |
516 490 | { |
517 491 | LogIO os(LogOrigin("ImagePolarimetry", "sigmaLinPolPosAng(...)", WHERE)); |
518 - | if (itsStokesPtr[ImagePolarimetry::Q]==0 && itsStokesPtr[ImagePolarimetry::U]==0) { |
492 + | if (_stokes[ImagePolarimetry::Q]==0 && _stokes[ImagePolarimetry::U]==0) { |
519 493 | os << "This image does not have Stokes Q and U so cannot provide linear polarization" << LogIO::EXCEPTION; |
520 494 | } |
521 495 | _checkQUBeams(false); |
522 496 | // Make expression |
523 497 | |
524 498 | Float sigma2 = 0.0; |
525 499 | if (sigma > 0) { |
526 500 | sigma2 = sigma; |
527 501 | } else { |
528 502 | sigma2 = ImagePolarimetry::sigma(clip); |
529 503 | } |
530 504 | Float fac = 0.5 * sigma2; |
531 505 | if (!radians) fac *= 180 / C::pi; |
532 506 | LatticeExprNode node(fac / |
533 - | amp(*itsStokesPtr[ImagePolarimetry::U], *itsStokesPtr[ImagePolarimetry::Q])); |
507 + | amp(*_stokes[ImagePolarimetry::U], *_stokes[ImagePolarimetry::Q])); |
534 508 | LatticeExpr<Float> le(node); |
535 509 | ImageExpr<Float> ie(le, String("LinearlyPolarizedPositionAngleError")); |
536 510 | // |
537 511 | if (radians) { |
538 512 | ie.setUnits(Unit("rad")); |
539 513 | } else { |
540 514 | ie.setUnits(Unit("deg")); |
541 515 | } |
542 - | ImageInfo ii = itsInImagePtr->imageInfo(); |
516 + | ImageInfo ii = _image->imageInfo(); |
543 517 | ii.removeRestoringBeam(); |
544 518 | ie.setImageInfo(ii); |
545 - | |
546 - | // Fiddle Stokes coordinate |
547 - | |
548 - | fiddleStokesCoordinate(ie, Stokes::Pangle); |
549 - | // |
519 + | _fiddleStokesCoordinate(ie, Stokes::Pangle); |
550 520 | return ie; |
551 521 | } |
552 522 | |
553 523 | Float ImagePolarimetry::sigma(Float clip) |
554 524 | { |
555 525 | LogIO os(LogOrigin("ImagePolarimetry", __func__, WHERE)); |
556 526 | Float sigma2 = 0.0; |
557 527 | |
558 - | if (itsStokesPtr[ImagePolarimetry::V]!=0) { |
528 + | if (_stokes[ImagePolarimetry::V]!=0) { |
559 529 | os << LogIO::NORMAL << "Determined noise from V image to be "; |
560 - | sigma2 = ImagePolarimetry::sigma(ImagePolarimetry::V, clip); |
530 + | sigma2 = _sigma(ImagePolarimetry::V, clip); |
561 531 | } |
562 532 | else if ( |
563 - | itsStokesPtr[ImagePolarimetry::Q]!=0 |
564 - | && itsStokesPtr[ImagePolarimetry::U]!=0 |
533 + | _stokes[ImagePolarimetry::Q]!=0 |
534 + | && _stokes[ImagePolarimetry::U]!=0 |
565 535 | && _checkQUBeams(false, false) |
566 536 | ) { |
567 537 | sigma2 = sigmaLinPolInt(clip); |
568 538 | } |
569 - | else if (itsStokesPtr[ImagePolarimetry::Q]!=0) { |
539 + | else if (_stokes[ImagePolarimetry::Q]!=0) { |
570 540 | os << LogIO::NORMAL << "Determined noise from Q image to be " << LogIO::POST; |
571 - | sigma2 = ImagePolarimetry::sigma(ImagePolarimetry::Q, clip); |
572 - | } else if (itsStokesPtr[ImagePolarimetry::U]!=0) { |
541 + | sigma2 = _sigma(ImagePolarimetry::Q, clip); |
542 + | } else if (_stokes[ImagePolarimetry::U]!=0) { |
573 543 | os << LogIO::NORMAL << "Determined noise from U image to be " << LogIO::POST; |
574 - | sigma2 = ImagePolarimetry::sigma(ImagePolarimetry::U, clip); |
575 - | } else if (itsStokesPtr[ImagePolarimetry::I]!=0) { |
544 + | sigma2 = _sigma(ImagePolarimetry::U, clip); |
545 + | } else if (_stokes[ImagePolarimetry::I]!=0) { |
576 546 | os << LogIO::NORMAL << "Determined noise from I image to be " << LogIO::POST; |
577 - | sigma2 = ImagePolarimetry::sigma(ImagePolarimetry::I, clip); |
547 + | sigma2 = _sigma(ImagePolarimetry::I, clip); |
578 548 | } |
579 549 | os << sigma2 << LogIO::POST; |
580 550 | return sigma2; |
581 551 | } |
582 552 | |
583 553 | void ImagePolarimetry::rotationMeasure( |
584 554 | ImageInterface<Float>*& rmOutPtr, ImageInterface<Float>*& rmOutErrorPtr, |
585 555 | ImageInterface<Float>*& pa0OutPtr, ImageInterface<Float>*& pa0OutErrorPtr, |
586 556 | ImageInterface<Float>*& nTurnsOutPtr, ImageInterface<Float>*& chiSqOutPtr, |
587 557 | Int axis, Float rmMax, Float maxPaErr, Float sigma, Float rmFg, |
588 558 | Bool showProgress |
589 559 | ) { |
590 560 | LogIO os(LogOrigin("ImagePolarimetry", __func__, WHERE)); |
591 - | hasQU(); |
561 + | _hasQU(); |
592 562 | _checkQUBeams(false); |
593 563 | // Do we have anything to do ? |
594 564 | if (!rmOutPtr && !rmOutErrorPtr && !pa0OutPtr && !pa0OutErrorPtr) { |
595 565 | os << "No output images specified" << LogIO::EXCEPTION; |
596 566 | } |
597 567 | // Find expected shape of output RM images (Stokes and spectral axes gone) |
598 568 | CoordinateSystem cSysRM; |
599 569 | Int fAxis, sAxis; |
600 570 | IPosition shapeRM = rotationMeasureShape(cSysRM, fAxis, sAxis, os, axis); |
601 571 | IPosition shapeNTurns = shapeRM; |
649 619 | } |
650 620 | // Do we have enough frequency pixels ? |
651 621 | const uInt nFreq = pa.shape()(fAxis); |
652 622 | if (nFreq < 3) { |
653 623 | os << "This image only has " << nFreq << "frequencies, this is not enough" |
654 624 | << LogIO::EXCEPTION; |
655 625 | } |
656 626 | // Copy units only over. The output images don't have a beam |
657 627 | // so unset beam. MiscInfo and history require writable II. |
658 628 | // We leave this to the caller who knows what sort of II these are. |
659 - | ImageInfo ii = itsInImagePtr->imageInfo(); |
629 + | ImageInfo ii = _image->imageInfo(); |
660 630 | ii.removeRestoringBeam(); |
661 631 | if (rmOutPtr) { |
662 632 | rmOutPtr->setImageInfo(ii); |
663 633 | rmOutPtr->setUnits(Unit("rad/m/m")); |
664 634 | } |
665 635 | if (rmOutErrorPtr) { |
666 636 | rmOutErrorPtr->setImageInfo(ii); |
667 637 | rmOutErrorPtr->setUnits(Unit("rad/m/m")); |
668 638 | } |
669 639 | if (pa0OutPtr) { |
699 669 | wsq(i) = csq / freqs(i) / freqs(i); // m**2 |
700 670 | } |
701 671 | // Sort into increasing wavelength |
702 672 | Vector<uInt> sortidx; |
703 673 | GenSortIndirect<Float>::sort( |
704 674 | sortidx, wsq, Sort::Ascending, Sort::QuickSort|Sort::NoDuplicates |
705 675 | ); |
706 676 | Vector<Float> wsqsort(sortidx.nelements()); |
707 677 | for (uInt i=0; i<wsqsort.nelements(); i++) wsqsort(i) = wsq(sortidx(i)); |
708 678 | // Make fitter |
709 - | if (itsFitterPtr==0) { |
710 - | itsFitterPtr = new LinearFitSVD<Float>; |
679 + | if (_fitter==0) { |
680 + | _fitter = new LinearFitSVD<Float>; |
711 681 | // Create and set the polynomial functional |
712 682 | // p = c(0) + c(1)*x where x = lambda**2 |
713 683 | // PA = PA0 + RM*Lambda**2 |
714 684 | Polynomial<AutoDiff<Float> > poly1(1); |
715 685 | // Makes a copy of poly1 |
716 - | itsFitterPtr->setFunction(poly1); |
686 + | _fitter->setFunction(poly1); |
717 687 | } |
718 688 | // Deal with masks. The outputs are all given a mask if possible as we |
719 689 | // don't know at this point whether output points will be masked or not |
720 690 | IPosition whereRM; |
721 691 | Bool isMaskedRM = false; |
722 692 | Lattice<Bool>* outRMMaskPtr = 0; |
723 693 | if (rmOutPtr) { |
724 - | isMaskedRM = dealWithMask (outRMMaskPtr, rmOutPtr, os, String("RM")); |
694 + | isMaskedRM = _dealWithMask (outRMMaskPtr, rmOutPtr, os, String("RM")); |
725 695 | whereRM.resize(rmOutPtr->ndim()); |
726 696 | whereRM = 0; |
727 697 | } |
728 698 | Bool isMaskedRMErr = false; |
729 699 | Lattice<Bool>* outRMErrMaskPtr = 0; |
730 700 | if (rmOutErrorPtr) { |
731 - | isMaskedRMErr = dealWithMask( |
701 + | isMaskedRMErr = _dealWithMask( |
732 702 | outRMErrMaskPtr, rmOutErrorPtr, os, String("RM error") |
733 703 | ); |
734 704 | whereRM.resize(rmOutErrorPtr->ndim()); |
735 705 | whereRM = 0; |
736 706 | } |
737 707 | IPosition wherePA; |
738 708 | Bool isMaskedPa0 = false; |
739 709 | Lattice<Bool>* outPa0MaskPtr = 0; |
740 710 | if (pa0OutPtr) { |
741 - | isMaskedPa0 = dealWithMask( |
711 + | isMaskedPa0 = _dealWithMask( |
742 712 | outPa0MaskPtr, pa0OutPtr, os, String("Position Angle") |
743 713 | ); |
744 714 | wherePA.resize(pa0OutPtr->ndim()); |
745 715 | wherePA = 0; |
746 716 | } |
747 717 | Bool isMaskedPa0Err = false; |
748 718 | Lattice<Bool>* outPa0ErrMaskPtr = 0; |
749 719 | if (pa0OutErrorPtr) { |
750 - | isMaskedPa0Err = dealWithMask( |
720 + | isMaskedPa0Err = _dealWithMask( |
751 721 | outPa0ErrMaskPtr, pa0OutErrorPtr, os, "Position Angle error" |
752 722 | ); |
753 723 | wherePA.resize(pa0OutErrorPtr->ndim()); |
754 724 | wherePA = 0; |
755 725 | } |
756 726 | IPosition whereNTurns; |
757 727 | Bool isMaskedNTurns = false; |
758 728 | Lattice<Bool>* outNTurnsMaskPtr = 0; |
759 729 | if (nTurnsOutPtr) { |
760 - | isMaskedNTurns = dealWithMask( |
730 + | isMaskedNTurns = _dealWithMask( |
761 731 | outNTurnsMaskPtr, nTurnsOutPtr, os, "nTurns" |
762 732 | ); |
763 733 | whereNTurns.resize(nTurnsOutPtr->ndim()); |
764 734 | whereNTurns = 0; |
765 735 | } |
766 736 | IPosition whereChiSq; |
767 737 | Bool isMaskedChiSq = false; |
768 738 | Lattice<Bool>* outChiSqMaskPtr = 0; |
769 739 | if (chiSqOutPtr) { |
770 - | isMaskedChiSq = dealWithMask( |
740 + | isMaskedChiSq = _dealWithMask( |
771 741 | outChiSqMaskPtr, chiSqOutPtr, os, String("chi sqared") |
772 742 | ); |
773 743 | whereChiSq.resize(chiSqOutPtr->ndim()); |
774 744 | whereChiSq = 0; |
775 745 | } |
776 746 | Array<Bool> tmpMaskRM(IPosition(shapeRM.nelements(), 1), true); |
777 747 | Array<Float> tmpValueRM(IPosition(shapeRM.nelements(), 1), 0.0); |
778 748 | Array<Bool> tmpMaskPA(IPosition(shapePA.nelements(), 1), true); |
779 749 | Array<Float> tmpValuePA(IPosition(shapePA.nelements(), 1), 0.0); |
780 750 | Array<Float> tmpValueNTurns(IPosition(shapeNTurns.nelements(), 1), 0.0); |
807 777 | nMin, nMax, "Profiles fitted", "Fitting", |
808 778 | "", "", true, max(1,Int(nMax/100)) |
809 779 | ) |
810 780 | ); |
811 781 | } |
812 782 | // As a (temporary?) workaround the cache of the main image is set up in |
813 783 | // such a way that it can hold the full frequency and stokes axes. |
814 784 | // The stokes axis is important, otherwise the cache is set up |
815 785 | // (by the TiledStMan) such that it can hold only 1 stokes |
816 786 | // with the result that iterating is tremendously slow. |
817 - | // We also need to cast the const away from itsInImagePtr. |
818 - | const IPosition mainShape = itsInImagePtr->shape(); |
787 + | // We also need to cast the const away from _image. |
788 + | const IPosition mainShape = _image->shape(); |
819 789 | uInt nrtiles = (1 + (mainShape(fAxis)-1) / tileShape(fAxis)) * |
820 790 | (1 + (mainShape(sAxis)-1) / tileShape(sAxis)); |
821 791 | ImageInterface<Float>* mainImagePtr = |
822 - | const_cast<ImageInterface<Float>*>(itsInImagePtr.get()); |
792 + | const_cast<ImageInterface<Float>*>(_image.get()); |
823 793 | mainImagePtr->setCacheSizeInTiles (nrtiles); |
824 794 | String posString; |
825 795 | Bool ok = false; |
826 796 | IPosition shp; |
827 797 | for (it.reset(); !it.atEnd(); it++) { |
828 798 | // Find rotation measure for this line |
829 - | ok = findRotationMeasure (rm, rmErr, pa0, pa0Err, rChiSq, nTurns, |
799 + | ok = _findRotationMeasure (rm, rmErr, pa0, pa0Err, rChiSq, nTurns, |
830 800 | sortidx, wsqsort, it.vectorCursor(), |
831 801 | it.getMask(false), |
832 802 | paerr.getSlice(it.position(),it.cursorShape()), |
833 803 | rmFg, rmMax, maxPaErr, /*plotter,*/ posString); |
834 804 | // Plonk values into output image. This is slow and clunky, but should be relatively fast |
835 805 | // c.f. the fitting. Could be reimplemented with LatticeApply if need be. Buffering |
836 806 | // is hard because the navigator doesn't take a regular path. If I used a LatticeStepper |
837 807 | // instead, the path would be regular and then I could buffer, but then the iteration |
838 808 | // would be less efficient !!! |
839 809 | |
919 889 | Int& sAxis, LogIO&, Int spectralAxis) const |
920 890 | { |
921 891 | |
922 892 | // Construction image CS |
923 893 | |
924 894 | CoordinateSystem cSys0 = coordinates(); |
925 895 | |
926 896 | // Find frequency axis |
927 897 | |
928 898 | Int spectralCoord; |
929 - | findFrequencyAxis (spectralCoord, fAxis, cSys0, spectralAxis); |
899 + | _findFrequencyAxis(spectralCoord, fAxis, cSys0, spectralAxis); |
930 900 | |
931 901 | // Find Stokes axis (we know it has one) |
932 902 | |
933 903 | Int afterCoord = -1; |
934 904 | Int stokesCoord = cSys0.findCoordinate(Coordinate::STOKES, afterCoord); |
935 905 | Vector<Int> pixelAxes = cSys0.pixelAxes(stokesCoord); |
936 906 | sAxis = pixelAxes(0); |
937 907 | |
938 908 | // What shape should the image be ? Frequency and stokes axes should be gone. |
939 909 | |
966 936 | Int& fAxis, Int& sAxis, LogIO&, Int spectralAxis) const |
967 937 | { |
968 938 | |
969 939 | // Construction image CS |
970 940 | |
971 941 | CoordinateSystem cSys0 = coordinates(); |
972 942 | |
973 943 | // Find frequency axis |
974 944 | |
975 945 | Int spectralCoord = -1; |
976 - | findFrequencyAxis (spectralCoord, fAxis, cSys0, spectralAxis); |
946 + | _findFrequencyAxis (spectralCoord, fAxis, cSys0, spectralAxis); |
977 947 | |
978 948 | // Find Stokes axis (we know it has one) |
979 949 | |
980 950 | Int afterCoord = -1; |
981 951 | Int stokesCoord = cSys0.findCoordinate(Coordinate::STOKES, afterCoord); |
982 952 | Vector<Int> pixelAxes = cSys0.pixelAxes(stokesCoord); |
983 953 | sAxis = pixelAxes(0); |
984 - | |
985 - | // Fiddle StokesCoordinate |
986 - | |
987 - | fiddleStokesCoordinate(cSys0, Stokes::Pangle); |
988 - | |
954 + | _fiddleStokesCoordinate(cSys0, Stokes::Pangle); |
989 955 | // Create output coordinate system |
990 956 | |
991 957 | CoordinateSystem tmp; |
992 958 | cSys = tmp; |
993 959 | for (Int i=0;i<Int(cSys0.nCoordinates()); i++) { |
994 960 | if (i!=spectralCoord) { |
995 961 | cSys.addCoordinate(cSys0.coordinate(i)); |
996 962 | } |
997 963 | } |
998 964 | |
1014 980 | } |
1015 981 | } |
1016 982 | } |
1017 983 | // |
1018 984 | return shape; |
1019 985 | } |
1020 986 | |
1021 987 | |
1022 988 | ImageExpr<Float> ImagePolarimetry::stokesI() const |
1023 989 | { |
1024 - | return makeStokesExpr(itsStokesPtr[ImagePolarimetry::I], Stokes::I, String("StokesI")); |
990 + | return _makeStokesExpr(_stokes[ImagePolarimetry::I], Stokes::I, String("StokesI")); |
1025 991 | } |
1026 992 | |
1027 993 | Float ImagePolarimetry::sigmaStokesI(Float clip) |
1028 994 | { |
1029 - | return ImagePolarimetry::sigma(ImagePolarimetry::I, clip); |
995 + | return _sigma(ImagePolarimetry::I, clip); |
1030 996 | } |
1031 997 | |
1032 998 | ImageExpr<Float> ImagePolarimetry::stokesQ() const |
1033 999 | { |
1034 - | return makeStokesExpr(itsStokesPtr[ImagePolarimetry::Q], Stokes::Q, String("StokesQ")); |
1000 + | return _makeStokesExpr(_stokes[ImagePolarimetry::Q], Stokes::Q, String("StokesQ")); |
1035 1001 | } |
1036 1002 | |
1037 1003 | Float ImagePolarimetry::sigmaStokesQ(Float clip) |
1038 1004 | { |
1039 - | return ImagePolarimetry::sigma(ImagePolarimetry::Q, clip); |
1005 + | return _sigma(ImagePolarimetry::Q, clip); |
1040 1006 | } |
1041 1007 | |
1042 1008 | ImageExpr<Float> ImagePolarimetry::stokesU() const |
1043 1009 | { |
1044 - | return makeStokesExpr(itsStokesPtr[ImagePolarimetry::U], Stokes::U, String("StokesU")); |
1010 + | return _makeStokesExpr(_stokes[ImagePolarimetry::U], Stokes::U, String("StokesU")); |
1045 1011 | } |
1046 1012 | |
1047 1013 | Float ImagePolarimetry::sigmaStokesU(Float clip) |
1048 1014 | { |
1049 - | return ImagePolarimetry::sigma(ImagePolarimetry::U, clip); |
1015 + | return _sigma(ImagePolarimetry::U, clip); |
1050 1016 | } |
1051 1017 | |
1052 1018 | ImageExpr<Float> ImagePolarimetry::stokesV() const |
1053 1019 | { |
1054 - | return makeStokesExpr(itsStokesPtr[ImagePolarimetry::V], Stokes::V, String("StokesV")); |
1020 + | return _makeStokesExpr(_stokes[ImagePolarimetry::V], Stokes::V, String("StokesV")); |
1055 1021 | } |
1056 1022 | |
1057 1023 | Float ImagePolarimetry::sigmaStokesV(Float clip) |
1058 1024 | { |
1059 - | return ImagePolarimetry::sigma(ImagePolarimetry::V, clip); |
1025 + | return _sigma(ImagePolarimetry::V, clip); |
1060 1026 | } |
1061 1027 | |
1062 1028 | ImageExpr<Float> ImagePolarimetry::stokes(ImagePolarimetry::StokesTypes stokes) const |
1063 1029 | { |
1064 - | Stokes::StokesTypes type = stokesType(stokes); |
1065 - | return makeStokesExpr(itsStokesPtr[stokes], type, stokesName(stokes)); |
1030 + | Stokes::StokesTypes type = _stokesType(stokes); |
1031 + | return _makeStokesExpr(_stokes[stokes], type, _stokesName(stokes)); |
1066 1032 | } |
1067 1033 | |
1068 1034 | Float ImagePolarimetry::sigmaStokes(ImagePolarimetry::StokesTypes stokes, Float clip) |
1069 1035 | { |
1070 - | return ImagePolarimetry::sigma(stokes, clip); |
1036 + | return _sigma(stokes, clip); |
1071 1037 | } |
1072 1038 | |
1073 1039 | void ImagePolarimetry::summary(LogIO& os) const |
1074 1040 | { |
1075 - | ImageSummary<Float> s(*itsInImagePtr); |
1041 + | ImageSummary<Float> s(*_image); |
1076 1042 | s.list(os); |
1077 1043 | } |
1078 1044 | |
1079 1045 | ImageExpr<Float> ImagePolarimetry::totPolInt(Bool debias, Float clip, Float sigma) |
1080 1046 | { |
1081 1047 | LogIO os(LogOrigin("ImagePolarimetry", __func__, WHERE)); |
1082 1048 | Bool doLin, doCirc; |
1083 1049 | _setDoLinDoCirc(doLin, doCirc); |
1084 1050 | // Make node. |
1085 1051 | |
1086 - | LatticeExprNode node = makePolIntNode(os, debias, clip, sigma, doLin, doCirc); |
1052 + | LatticeExprNode node = _makePolIntNode(os, debias, clip, sigma, doLin, doCirc); |
1087 1053 | |
1088 1054 | // Make expression |
1089 1055 | |
1090 1056 | LatticeExpr<Float> le(node); |
1091 1057 | ImageExpr<Float> ie(le, String("totalPolarizedIntensity")); |
1092 - | ie.setUnits(itsInImagePtr->units()); // Dodgy. The beam is now rectified |
1093 - | StokesTypes stokes = itsStokesPtr[Q] != 0 |
1058 + | ie.setUnits(_image->units()); // Dodgy. The beam is now rectified |
1059 + | StokesTypes stokes = _stokes[Q] != 0 |
1094 1060 | ? Q |
1095 - | : itsStokesPtr[U] != 0 |
1061 + | : _stokes[U] != 0 |
1096 1062 | ? U |
1097 1063 | : V; |
1098 1064 | _setInfo(ie, stokes); |
1099 - | |
1100 - | // Fiddle Stokes coordinate in ImageExpr |
1101 - | |
1102 - | fiddleStokesCoordinate(ie, Stokes::Ptotal); |
1103 - | // |
1065 + | _fiddleStokesCoordinate(ie, Stokes::Ptotal); |
1104 1066 | return ie; |
1105 1067 | } |
1106 1068 | |
1107 1069 | Float ImagePolarimetry::sigmaTotPolInt(Float clip, Float sigma) |
1108 1070 | // |
1109 1071 | // sigma_P = sigma_QUV |
1110 1072 | // |
1111 1073 | { |
1112 1074 | LogIO os(LogOrigin("ImagePolarimetry", __func__, WHERE)); |
1113 1075 | Bool doLin, doCirc; |
1123 1085 | } |
1124 1086 | return sigma2; |
1125 1087 | } |
1126 1088 | |
1127 1089 | |
1128 1090 | IPosition ImagePolarimetry::singleStokesShape(CoordinateSystem& cSys, Stokes::StokesTypes type) const |
1129 1091 | { |
1130 1092 | // We know the image has a Stokes coordinate or it |
1131 1093 | // would have failed at construction |
1132 1094 | |
1133 - | CoordinateSystem cSys0 = itsInImagePtr->coordinates(); |
1134 - | fiddleStokesCoordinate(cSys0, type); |
1095 + | CoordinateSystem cSys0 = _image->coordinates(); |
1096 + | _fiddleStokesCoordinate(cSys0, type); |
1135 1097 | cSys = cSys0; |
1136 1098 | // |
1137 1099 | Int afterCoord = -1; |
1138 1100 | Int iStokes = cSys0.findCoordinate(Coordinate::STOKES, afterCoord); |
1139 1101 | Vector<Int> pixelAxes = cSys0.pixelAxes(iStokes); |
1140 - | IPosition shape = itsInImagePtr->shape(); |
1102 + | IPosition shape = _image->shape(); |
1141 1103 | shape(pixelAxes(0)) = 1; |
1142 1104 | // |
1143 1105 | return shape; |
1144 1106 | } |
1145 1107 | |
1146 1108 | |
1147 1109 | ImageExpr<Float> ImagePolarimetry::depolarizationRatio (const ImageInterface<Float>& im1, |
1148 1110 | const ImageInterface<Float>& im2, |
1149 1111 | Bool debias, Float clip, Float sigma) |
1150 1112 | { |
1182 1144 | LatticeExpr<Float> le(n3); |
1183 1145 | ImageExpr<Float> sigmaDepol(le, String("DepolarizationRatioError")); |
1184 1146 | return sigmaDepol; |
1185 1147 | } |
1186 1148 | |
1187 1149 | |
1188 1150 | |
1189 1151 | // Private functions |
1190 1152 | |
1191 1153 | |
1192 - | void ImagePolarimetry::cleanup() |
1154 + | void ImagePolarimetry::_cleanup() |
1193 1155 | { |
1194 - | itsInImagePtr.reset(); |
1156 + | _image.reset(); |
1195 1157 | // |
1196 1158 | for (uInt i=0; i<4; i++) { |
1197 - | delete itsStokesPtr[i]; |
1198 - | itsStokesPtr[i] = 0; |
1159 + | delete _stokes[i]; |
1160 + | _stokes[i] = 0; |
1199 1161 | // |
1200 - | delete itsStokesStatsPtr[i]; |
1201 - | itsStokesStatsPtr[i] = 0; |
1162 + | delete _stokesStats[i]; |
1163 + | _stokesStats[i] = 0; |
1202 1164 | } |
1203 1165 | // |
1204 - | if (itsFitterPtr!= 0) { |
1205 - | delete itsFitterPtr; |
1206 - | itsFitterPtr = 0; |
1166 + | if (_fitter!= 0) { |
1167 + | delete _fitter; |
1168 + | _fitter = 0; |
1207 1169 | } |
1208 1170 | } |
1209 1171 | |
1210 1172 | |
1211 - | void ImagePolarimetry::copyDataAndMask(ImageInterface<Float>& out, |
1173 + | void ImagePolarimetry::_copyDataAndMask(ImageInterface<Float>& out, |
1212 1174 | ImageInterface<Float>& in) const |
1213 1175 | // |
1214 1176 | // Copy the data and mask from input to output. If the input is |
1215 1177 | // masked, the output must already be masked and ready |
1216 1178 | // |
1217 1179 | { |
1218 1180 | // Do we need to stuff about with masks ? |
1219 1181 | |
1220 1182 | Bool doMask = in.isMasked() && out.hasPixelMask(); |
1221 1183 | Lattice<Bool>* pMaskOut = 0; |
1238 1200 | RO_MaskedLatticeIterator<Float> iter(in, stepper); |
1239 1201 | for (iter.reset(); !iter.atEnd(); iter++) { |
1240 1202 | out.putSlice (iter.cursor(), iter.position()); |
1241 1203 | if (doMask) { |
1242 1204 | pMaskOut->putSlice(iter.getMask(false), iter.position()); |
1243 1205 | } |
1244 1206 | } |
1245 1207 | } |
1246 1208 | |
1247 1209 | |
1248 - | void ImagePolarimetry::findFrequencyAxis (Int& spectralCoord, Int& fAxis, |
1210 + | void ImagePolarimetry::_findFrequencyAxis (Int& spectralCoord, Int& fAxis, |
1249 1211 | const CoordinateSystem& cSys, Int spectralAxis) const |
1250 1212 | { |
1251 - | LogIO os(LogOrigin("ImagePolarimetry", "findFrequencyAxis(...)", WHERE)); |
1213 + | LogIO os(LogOrigin("ImagePolarimetry", __func__, WHERE)); |
1252 1214 | spectralCoord = -1; |
1253 1215 | fAxis = -1; |
1254 1216 | if (spectralAxis >=0) { |
1255 1217 | if (spectralAxis < Int(cSys.nPixelAxes())) { |
1256 1218 | fAxis = spectralAxis; |
1257 1219 | Int axisInCoordinate; |
1258 1220 | cSys.findPixelAxis(spectralCoord, axisInCoordinate, fAxis); |
1259 1221 | |
1260 1222 | // Check coordinate type is one of expected types |
1261 1223 | |
1262 1224 | Bool ok = cSys.type(spectralCoord)==Coordinate::TABULAR || |
1263 1225 | cSys.type(spectralCoord)==Coordinate::LINEAR || |
1264 1226 | cSys.type(spectralCoord)==Coordinate::SPECTRAL; |
1265 1227 | if (!ok) { |
1266 1228 | os << "The specified axis of type " << cSys.showType(spectralCoord) |
1267 1229 | << " cannot be a frequency axis" << LogIO::EXCEPTION; |
1268 1230 | } |
1269 1231 | } else { |
1270 1232 | os << "Illegal spectral axis " << spectralAxis+1 << " given" << LogIO::EXCEPTION; |
1271 1233 | } |
1272 1234 | } else { |
1273 - | spectralCoord = findSpectralCoordinate(cSys, os, false); |
1235 + | spectralCoord = _findSpectralCoordinate(cSys, os, false); |
1274 1236 | if (spectralCoord < 0) { |
1275 1237 | for (uInt i=0; i<cSys.nCoordinates(); i++) { |
1276 1238 | if (cSys.type(i)==Coordinate::TABULAR || |
1277 1239 | cSys.type(i)==Coordinate::LINEAR) { |
1278 1240 | Vector<String> axisNames = cSys.coordinate(i).worldAxisNames(); |
1279 1241 | String tmp = axisNames(0); |
1280 1242 | tmp.upcase(); |
1281 1243 | if (tmp.contains(String("FREQ"))) { |
1282 1244 | spectralCoord = i; |
1283 1245 | break; |
1288 1250 | if (spectralCoord < 0) { |
1289 1251 | os << "Cannot find SpectralCoordinate in this image" << LogIO::EXCEPTION; |
1290 1252 | } else { |
1291 1253 | Vector<Int> pixelAxes = cSys.pixelAxes(spectralCoord); |
1292 1254 | fAxis = pixelAxes(0); |
1293 1255 | } |
1294 1256 | } |
1295 1257 | } |
1296 1258 | |
1297 1259 | |
1298 - | void ImagePolarimetry::findStokes() |
1260 + | void ImagePolarimetry::_findStokes() |
1299 1261 | { |
1300 - | LogIO os(LogOrigin("ImagePolarimetry", "findStokes(...)", WHERE)); |
1262 + | LogIO os(LogOrigin("ImagePolarimetry", __func__, WHERE)); |
1301 1263 | |
1302 1264 | // Do we have any Stokes ? |
1303 1265 | |
1304 - | const CoordinateSystem& cSys = itsInImagePtr->coordinates(); |
1266 + | const CoordinateSystem& cSys = _image->coordinates(); |
1305 1267 | Int afterCoord = -1; |
1306 1268 | Int iStokes = cSys.findCoordinate(Coordinate::STOKES, afterCoord); |
1307 1269 | if (iStokes<0) { |
1308 - | cleanup(); |
1270 + | _cleanup(); |
1309 1271 | os << "There is no Stokes Coordinate in this image" << LogIO::EXCEPTION; |
1310 1272 | } |
1311 1273 | if (afterCoord>0) { |
1312 1274 | os << LogIO::WARN |
1313 1275 | << "There is more than one Stokes coordinate in this image. Only first considered" << LogIO::POST; |
1314 1276 | } |
1315 1277 | |
1316 1278 | // Find the pixel axis of the image which is Stokes |
1317 1279 | |
1318 1280 | Vector<Int> pixelAxes = cSys.pixelAxes(iStokes); |
1319 1281 | |
1320 1282 | // Make the regions |
1321 1283 | |
1322 1284 | const StokesCoordinate& stokes = cSys.stokesCoordinate(iStokes); |
1323 - | const uInt ndim = itsInImagePtr->ndim(); |
1324 - | IPosition shape = itsInImagePtr->shape(); |
1285 + | const uInt ndim = _image->ndim(); |
1286 + | IPosition shape = _image->shape(); |
1325 1287 | IPosition blc(ndim,0); |
1326 1288 | IPosition trc(shape-1); |
1327 1289 | // |
1328 1290 | Int pix; |
1329 1291 | if (stokes.toPixel(pix, Stokes::I)) { |
1330 - | itsStokesPtr[ImagePolarimetry::I] = makeSubImage(blc, trc, pixelAxes(0), pix); |
1292 + | _stokes[ImagePolarimetry::I] = _makeSubImage(blc, trc, pixelAxes(0), pix); |
1331 1293 | } |
1332 1294 | if (stokes.toPixel(pix, Stokes::Q)) { |
1333 - | itsStokesPtr[ImagePolarimetry::Q] = makeSubImage(blc, trc, pixelAxes(0), pix); |
1295 + | _stokes[ImagePolarimetry::Q] = _makeSubImage(blc, trc, pixelAxes(0), pix); |
1334 1296 | } |
1335 1297 | if (stokes.toPixel(pix, Stokes::U)) { |
1336 - | itsStokesPtr[ImagePolarimetry::U] = makeSubImage(blc, trc, pixelAxes(0), pix); |
1298 + | _stokes[ImagePolarimetry::U] = _makeSubImage(blc, trc, pixelAxes(0), pix); |
1337 1299 | } |
1338 1300 | if (stokes.toPixel(pix, Stokes::V)) { |
1339 - | itsStokesPtr[ImagePolarimetry::V] = makeSubImage(blc, trc, pixelAxes(0), pix); |
1301 + | _stokes[ImagePolarimetry::V] = _makeSubImage(blc, trc, pixelAxes(0), pix); |
1340 1302 | } |
1341 1303 | // |
1342 - | if ( (itsStokesPtr[ImagePolarimetry::Q]!=0 && itsStokesPtr[ImagePolarimetry::U]==0) || |
1343 - | (itsStokesPtr[ImagePolarimetry::Q]==0 && itsStokesPtr[ImagePolarimetry::U]!=0)) { |
1344 - | cleanup(); |
1304 + | if ( (_stokes[ImagePolarimetry::Q]!=0 && _stokes[ImagePolarimetry::U]==0) || |
1305 + | (_stokes[ImagePolarimetry::Q]==0 && _stokes[ImagePolarimetry::U]!=0)) { |
1306 + | _cleanup(); |
1345 1307 | os << "This Stokes coordinate has only one of Q and U. This is not useful" << LogIO::EXCEPTION; |
1346 1308 | } |
1347 - | if (itsStokesPtr[ImagePolarimetry::Q]==0 && |
1348 - | itsStokesPtr[ImagePolarimetry::U]==0 && |
1349 - | itsStokesPtr[ImagePolarimetry::V]==0) { |
1350 - | cleanup(); |
1309 + | if (_stokes[ImagePolarimetry::Q]==0 && |
1310 + | _stokes[ImagePolarimetry::U]==0 && |
1311 + | _stokes[ImagePolarimetry::V]==0) { |
1312 + | _cleanup(); |
1351 1313 | os << "This image has no Stokes Q, U, or V. This is not useful" << LogIO::EXCEPTION; |
1352 1314 | } |
1353 1315 | } |
1354 1316 | |
1355 1317 | |
1356 - | void ImagePolarimetry::fiddleStokesCoordinate(ImageInterface<Float>& im, Stokes::StokesTypes type) const |
1318 + | void ImagePolarimetry::_fiddleStokesCoordinate(ImageInterface<Float>& im, Stokes::StokesTypes type) const |
1357 1319 | { |
1358 1320 | CoordinateSystem cSys = im.coordinates(); |
1359 - | fiddleStokesCoordinate(cSys, type); |
1321 + | _fiddleStokesCoordinate(cSys, type); |
1360 1322 | im.setCoordinateInfo(cSys); |
1361 1323 | } |
1362 1324 | |
1363 - | void ImagePolarimetry::fiddleStokesCoordinate(CoordinateSystem& cSys, Stokes::StokesTypes type) const |
1325 + | void ImagePolarimetry::_fiddleStokesCoordinate(CoordinateSystem& cSys, Stokes::StokesTypes type) const |
1364 1326 | { |
1365 1327 | Int afterCoord = -1; |
1366 1328 | Int iStokes = cSys.findCoordinate(Coordinate::STOKES, afterCoord); |
1367 1329 | // |
1368 1330 | Vector<Int> which(1); |
1369 1331 | which(0) = Int(type); |
1370 1332 | StokesCoordinate stokes(which); |
1371 1333 | cSys.replaceCoordinate(stokes, iStokes); |
1372 1334 | } |
1373 1335 | |
1374 - | void ImagePolarimetry::fiddleStokesCoordinate(ImageInterface<Complex>& ie, Stokes::StokesTypes type) const |
1336 + | void ImagePolarimetry::_fiddleStokesCoordinate(ImageInterface<Complex>& ie, Stokes::StokesTypes type) const |
1375 1337 | { |
1376 1338 | CoordinateSystem cSys = ie.coordinates(); |
1377 1339 | // |
1378 1340 | Int afterCoord = -1; |
1379 1341 | Int iStokes = cSys.findCoordinate(Coordinate::STOKES, afterCoord); |
1380 1342 | // |
1381 1343 | Vector<Int> which(1); |
1382 1344 | which(0) = Int(type); |
1383 1345 | StokesCoordinate stokes(which); |
1384 1346 | cSys.replaceCoordinate(stokes, iStokes); |
1385 1347 | ie.setCoordinateInfo(cSys); |
1386 1348 | } |
1387 1349 | |
1388 - | void ImagePolarimetry::fiddleTimeCoordinate(ImageInterface<Complex>& ie, const Quantum<Double>& f, |
1350 + | void ImagePolarimetry::_fiddleTimeCoordinate(ImageInterface<Complex>& ie, const Quantum<Double>& f, |
1389 1351 | Int coord) const |
1390 1352 | { |
1391 - | LogIO os(LogOrigin("ImagePolarimetry", "fiddleTimeCoordinate(...)", WHERE)); |
1353 + | LogIO os(LogOrigin("ImagePolarimetry", __func__, WHERE)); |
1392 1354 | // |
1393 1355 | CoordinateSystem cSys = ie.coordinates(); |
1394 1356 | Coordinate* pC = cSys.coordinate(coord).clone(); |
1395 1357 | AlwaysAssert(pC->nPixelAxes()==1,AipsError); |
1396 1358 | AlwaysAssert(pC->type()==Coordinate::LINEAR,AipsError); |
1397 1359 | // |
1398 1360 | Vector<String> axisUnits = pC->worldAxisUnits(); |
1399 1361 | axisUnits = String("s"); |
1400 1362 | if (!pC->setWorldAxisUnits(axisUnits)) { |
1401 1363 | os << "Failed to set TimeCoordinate units to seconds because " << pC->errorMessage() << LogIO::EXCEPTION; |
1416 1378 | // |
1417 1379 | LinearCoordinate lC(axisNames, axisUnits, refVal, inc, |
1418 1380 | pC->linearTransform().copy(), pC->referencePixel().copy()); |
1419 1381 | // |
1420 1382 | cSys.replaceCoordinate(lC, coord); |
1421 1383 | ie.setCoordinateInfo(cSys); |
1422 1384 | delete pC; |
1423 1385 | } |
1424 1386 | |
1425 1387 | |
1426 - | Quantum<Double> ImagePolarimetry::findCentralFrequency(const Coordinate& coord, Int shape) const |
1388 + | Quantum<Double> ImagePolarimetry::_findCentralFrequency(const Coordinate& coord, Int shape) const |
1427 1389 | { |
1428 1390 | AlwaysAssert(coord.nPixelAxes()==1,AipsError); |
1429 1391 | // |
1430 1392 | Vector<Double> pixel(1); |
1431 1393 | Vector<Double> world; |
1432 1394 | pixel(0) = Double(shape - 1) / 2.0; |
1433 1395 | if (!coord.toWorld(world, pixel)) { |
1434 - | LogIO os(LogOrigin("ImagePolarimetry", "findCentralFrequency(...)", WHERE)); |
1396 + | LogIO os(LogOrigin("ImagePolarimetry", __func__, WHERE)); |
1435 1397 | os << "Failed to convert pixel to world for SpectralCoordinate because " |
1436 1398 | << coord.errorMessage() << LogIO::EXCEPTION; |
1437 1399 | } |
1438 1400 | Vector<String> units = coord.worldAxisUnits(); |
1439 1401 | return Quantum<Double>(world(0), units(0)); |
1440 1402 | } |
1441 1403 | |
1442 1404 | |
1443 - | Int ImagePolarimetry::findSpectralCoordinate(const CoordinateSystem& cSys, LogIO& os, |
1405 + | Int ImagePolarimetry::_findSpectralCoordinate(const CoordinateSystem& cSys, LogIO& os, |
1444 1406 | Bool fail) const |
1445 1407 | { |
1446 1408 | Int afterCoord = -1; |
1447 1409 | Int coord = cSys.findCoordinate(Coordinate::SPECTRAL, afterCoord); |
1448 1410 | if (coord<0) { |
1449 1411 | if (fail) os << "No spectral coordinate in this image" << LogIO::EXCEPTION; |
1450 1412 | } |
1451 1413 | if (afterCoord>0) { |
1452 1414 | os << LogIO::WARN << "This image has more than one spectral coordinate; only first used" |
1453 1415 | << LogIO::POST; |
1454 1416 | } |
1455 1417 | return coord; |
1456 1418 | } |
1457 1419 | |
1458 - | Bool ImagePolarimetry::findRotationMeasure (Float& rmFitted, Float& rmErrFitted, |
1420 + | Bool ImagePolarimetry::_findRotationMeasure (Float& rmFitted, Float& rmErrFitted, |
1459 1421 | Float& pa0Fitted, Float& pa0ErrFitted, |
1460 1422 | Float& rChiSqFitted, Float& nTurns, |
1461 1423 | const Vector<uInt>& sortidx, |
1462 1424 | const Vector<Float>& wsq2, const Vector<Float>& pa2, |
1463 1425 | const Array<Bool>& paMask2, |
1464 1426 | const Array<Float>& paerr2, |
1465 1427 | Float rmFg, Float rmMax, Float maxPaErr, |
1466 1428 | const String& posString) |
1467 1429 | // |
1468 1430 | // wsq is lambda squared in m**2 in increasing wavelength order |
1505 1467 | if (n<=1) return false; |
1506 1468 | // |
1507 1469 | pa.resize(n,true); |
1508 1470 | paerr.resize(n,true); |
1509 1471 | wsq.resize(n, true); |
1510 1472 | |
1511 1473 | // Treat supplementary and primary points separately |
1512 1474 | |
1513 1475 | Bool ok = false; |
1514 1476 | if (n==2) { |
1515 - | ok = rmSupplementaryFit(nTurns, rmFitted, rmErrFitted, pa0Fitted, pa0ErrFitted, |
1477 + | ok = _rmSupplementaryFit(nTurns, rmFitted, rmErrFitted, pa0Fitted, pa0ErrFitted, |
1516 1478 | rChiSqFitted, wsq, pa, paerr); |
1517 1479 | } else { |
1518 - | ok = rmPrimaryFit(nTurns, rmFitted, rmErrFitted, pa0Fitted, pa0ErrFitted, |
1480 + | ok = _rmPrimaryFit(nTurns, rmFitted, rmErrFitted, pa0Fitted, pa0ErrFitted, |
1519 1481 | rChiSqFitted, wsq, pa, paerr, rmMax, /*plotter,*/ posString); |
1520 1482 | } |
1521 1483 | |
1522 1484 | // Put position angle into the range 0->pi |
1523 1485 | |
1524 1486 | static MVAngle tmpMVA1; |
1525 1487 | if (ok) { |
1526 1488 | MVAngle tmpMVA0(pa0Fitted); |
1527 1489 | tmpMVA1 = tmpMVA0.binorm(0.0); |
1528 1490 | pa0Fitted = tmpMVA1.radian(); |
1529 1491 | |
1530 1492 | // Add foreground back on |
1531 1493 | |
1532 1494 | rmFitted += rmFg; |
1533 1495 | } |
1534 1496 | return ok; |
1535 1497 | } |
1536 1498 | |
1537 - | void ImagePolarimetry::hasQU () const |
1499 + | void ImagePolarimetry::_hasQU () const |
1538 1500 | { |
1539 - | Bool has = itsStokesPtr[ImagePolarimetry::Q]!=0 && |
1540 - | itsStokesPtr[ImagePolarimetry::U]!=0; |
1501 + | Bool has = _stokes[ImagePolarimetry::Q]!=0 && |
1502 + | _stokes[ImagePolarimetry::U]!=0; |
1541 1503 | if (!has) { |
1542 - | LogIO os(LogOrigin("ImagePolarimetry", "hasQU(...)", WHERE)); |
1504 + | LogIO os(LogOrigin("ImagePolarimetry", __func__, WHERE)); |
1543 1505 | os << "This image does not have Stokes Q and U which are required for this function" << LogIO::EXCEPTION; |
1544 1506 | } |
1545 1507 | } |
1546 1508 | |
1547 1509 | |
1548 - | ImageExpr<Float> ImagePolarimetry::makeStokesExpr(ImageInterface<Float>* imPtr, |
1510 + | ImageExpr<Float> ImagePolarimetry::_makeStokesExpr(ImageInterface<Float>* imPtr, |
1549 1511 | Stokes::StokesTypes type, const String& name) const |
1550 1512 | { |
1551 - | LogIO os(LogOrigin("ImagePolarimetry", "makeStokesExpr(...)", WHERE)); |
1513 + | LogIO os(LogOrigin("ImagePolarimetry", __func__, WHERE)); |
1552 1514 | if (imPtr==0) { |
1553 1515 | os << "This image does not have Stokes " << Stokes::name(type) << LogIO::EXCEPTION; |
1554 1516 | } |
1555 1517 | |
1556 1518 | // Make node. |
1557 1519 | |
1558 1520 | LatticeExprNode node(*imPtr); |
1559 1521 | |
1560 1522 | // Make expression |
1561 1523 | |
1562 1524 | LatticeExpr<Float> le(node); |
1563 1525 | ImageExpr<Float> ie(le, name); |
1564 - | ie.setUnits(itsInImagePtr->units()); |
1565 - | fiddleStokesCoordinate(ie, type); |
1566 - | // |
1526 + | ie.setUnits(_image->units()); |
1527 + | _fiddleStokesCoordinate(ie, type); |
1528 + | |
1567 1529 | return ie; |
1568 1530 | } |
1569 1531 | |
1570 1532 | |
1571 1533 | |
1572 - | ImageInterface<Float>* ImagePolarimetry::makeSubImage (IPosition& blc, |
1534 + | ImageInterface<Float>* ImagePolarimetry::_makeSubImage (IPosition& blc, |
1573 1535 | IPosition& trc, |
1574 1536 | Int axis, Int pix) const |
1575 1537 | { |
1576 1538 | blc(axis) = pix; |
1577 1539 | trc(axis) = pix; |
1578 1540 | LCSlicer slicer(blc, trc, RegionType::Abs); |
1579 1541 | ImageRegion region(slicer); |
1580 - | return new SubImage<Float>(*itsInImagePtr, region); |
1542 + | return new SubImage<Float>(*_image, region); |
1581 1543 | } |
1582 1544 | |
1583 1545 | |
1584 - | LatticeExprNode ImagePolarimetry::makePolIntNode(LogIO& os, Bool debias, Float clip, Float sigma, |
1546 + | LatticeExprNode ImagePolarimetry::_makePolIntNode(LogIO& os, Bool debias, Float clip, Float sigma, |
1585 1547 | Bool doLin, Bool doCirc) |
1586 1548 | { |
1587 1549 | LatticeExprNode linNode, circNode, node; |
1588 1550 | // |
1589 1551 | Float sigma2 = 0.0; |
1590 1552 | if (doLin) { |
1591 1553 | if (debias) { |
1592 1554 | if (sigma > 0.0) { |
1593 1555 | sigma2 = sigma; |
1594 1556 | } else { |
1595 1557 | sigma2 = ImagePolarimetry::sigma(clip); |
1596 1558 | } |
1597 1559 | } |
1598 - | linNode = LatticeExprNode(pow(*itsStokesPtr[ImagePolarimetry::U],2) + |
1599 - | pow(*itsStokesPtr[ImagePolarimetry::Q],2)); |
1560 + | linNode = LatticeExprNode(pow(*_stokes[ImagePolarimetry::U],2) + |
1561 + | pow(*_stokes[ImagePolarimetry::Q],2)); |
1600 1562 | } |
1601 1563 | // |
1602 1564 | if (doCirc) { |
1603 1565 | if (debias) { |
1604 1566 | if (sigma > 0.0) { |
1605 1567 | sigma2 = sigma; |
1606 1568 | } else { |
1607 1569 | sigma2 = ImagePolarimetry::sigma(clip); |
1608 1570 | } |
1609 1571 | } |
1610 - | circNode = LatticeExprNode(pow(*itsStokesPtr[ImagePolarimetry::V],2)); |
1572 + | circNode = LatticeExprNode(pow(*_stokes[ImagePolarimetry::V],2)); |
1611 1573 | } |
1612 1574 | // |
1613 1575 | Float sigmasq = sigma2 * sigma2; |
1614 1576 | if (doLin && doCirc) { |
1615 1577 | if (debias) { |
1616 1578 | node = linNode + circNode - LatticeExprNode(sigmasq); |
1617 1579 | os << LogIO::NORMAL << "Debiasing with sigma = " << sqrt(sigmasq) << LogIO::POST; |
1618 1580 | } else { |
1619 1581 | node = linNode + circNode; |
1620 1582 | } |
1631 1593 | os << LogIO::NORMAL << "Debiasing with sigma = " << sqrt(sigmasq) << LogIO::POST; |
1632 1594 | } else { |
1633 1595 | node = circNode; |
1634 1596 | } |
1635 1597 | } |
1636 1598 | // |
1637 1599 | return LatticeExprNode(sqrt(node)); |
1638 1600 | } |
1639 1601 | |
1640 1602 | |
1641 - | Bool ImagePolarimetry::rmPrimaryFit(Float& nTurns, Float& rmFitted, Float& rmErrFitted, |
1603 + | Bool ImagePolarimetry::_rmPrimaryFit(Float& nTurns, Float& rmFitted, Float& rmErrFitted, |
1642 1604 | Float& pa0Fitted, Float& pa0ErrFitted, |
1643 1605 | Float& rChiSqFitted, const Vector<Float>& wsq, |
1644 1606 | const Vector<Float>& pa, const Vector<Float>& paerr, |
1645 1607 | Float rmMax, const String& /*posString*/) |
1646 1608 | { |
1647 1609 | static Vector<Float> plotPA; |
1648 1610 | static Vector<Float> plotPAErr; |
1649 1611 | static Vector<Float> plotPAErrY1; |
1650 1612 | static Vector<Float> plotPAErrY2; |
1651 1613 | static Vector<Float> plotPAFit; |
1696 1658 | // |
1697 1659 | t = 0.5; |
1698 1660 | if (diff < 0) t = -0.5; |
1699 1661 | Int npi = Int(diff/C::pi + t); |
1700 1662 | fitpa(k) = pa(k) + npi*C::pi; |
1701 1663 | } |
1702 1664 | fitpa(0) = pa(0); |
1703 1665 | |
1704 1666 | // Do least squares fit |
1705 1667 | |
1706 - | if (!rmLsqFit (pars, wsq, fitpa, paerr)) return false; |
1668 + | if (!_rmLsqFit (pars, wsq, fitpa, paerr)) return false; |
1707 1669 | |
1708 1670 | // |
1709 1671 | if (pars(4) < chiSq) { |
1710 1672 | chiSq = pars(4); |
1711 1673 | // |
1712 1674 | nTurns = h; // Number of turns |
1713 1675 | rmFitted = pars(0); // Fitted RM |
1714 1676 | rmErrFitted = pars(1); // Error in RM |
1715 1677 | pa0Fitted = pars(2); // Fitted intrinsic angle |
1716 1678 | pa0ErrFitted = pars(3); // Error in angle |
1761 1723 | plotter.line(wsq, plotPAFit); |
1762 1724 | |
1763 1725 | } |
1764 1726 | */ |
1765 1727 | // |
1766 1728 | return true; |
1767 1729 | } |
1768 1730 | |
1769 1731 | |
1770 1732 | |
1771 - | Bool ImagePolarimetry::rmSupplementaryFit(Float& nTurns, Float& rmFitted, Float& rmErrFitted, |
1733 + | Bool ImagePolarimetry::_rmSupplementaryFit(Float& nTurns, Float& rmFitted, Float& rmErrFitted, |
1772 1734 | Float& pa0Fitted, Float& pa0ErrFitted, |
1773 1735 | Float& rChiSqFitted, const Vector<Float>& wsq, |
1774 1736 | const Vector<Float>& pa, const Vector<Float>& paerr) |
1775 1737 | { |
1776 1738 | |
1777 1739 | // For supplementary points find lowest residual RM |
1778 1740 | |
1779 1741 | const uInt n = wsq.nelements(); |
1780 1742 | // |
1781 1743 | Float absRM = 1e30; |
1782 1744 | Vector<Float> fitpa(pa.copy()); |
1783 1745 | Vector<Float> pars; |
1784 1746 | for (Int i=-2; i<3; i++) { |
1785 1747 | fitpa(n-1) = pa(n-1) + C::pi*i; |
1786 1748 | |
1787 1749 | // Do least squares fit |
1788 1750 | |
1789 - | if (!rmLsqFit (pars, wsq, fitpa, paerr)) return false; |
1751 + | if (! _rmLsqFit (pars, wsq, fitpa, paerr)) return false; |
1790 1752 | |
1791 1753 | // Save solution with lowest absolute RM |
1792 1754 | |
1793 1755 | if (abs(pars(0)) < absRM) { |
1794 1756 | absRM = abs(pars(0)); |
1795 1757 | // |
1796 1758 | nTurns = i; // nTurns |
1797 1759 | rmFitted = pars(0); // Fitted RM |
1798 1760 | rmErrFitted = pars(1); // Error in RM |
1799 1761 | pa0Fitted = pars(2); // Fitted intrinsic angle |
1803 1765 | } |
1804 1766 | |
1805 1767 | } |
1806 1768 | // |
1807 1769 | return true; |
1808 1770 | } |
1809 1771 | |
1810 1772 | |
1811 1773 | |
1812 1774 | |
1813 - | Bool ImagePolarimetry::rmLsqFit (Vector<Float>& pars, const Vector<Float>& wsq, |
1775 + | Bool ImagePolarimetry::_rmLsqFit (Vector<Float>& pars, const Vector<Float>& wsq, |
1814 1776 | const Vector<Float> pa, const Vector<Float>& paerr) const |
1815 1777 | { |
1816 1778 | |
1817 1779 | // Perform fit on unmasked data |
1818 1780 | |
1819 1781 | static Vector<Float> solution; |
1820 1782 | try { |
1821 - | solution = itsFitterPtr->fit(wsq, pa, paerr); |
1783 + | solution = _fitter->fit(wsq, pa, paerr); |
1822 1784 | } catch (AipsError x) { |
1823 1785 | return false; |
1824 1786 | } |
1825 1787 | // |
1826 - | const Vector<Double>& cv = itsFitterPtr->compuCovariance().diagonal(); |
1788 + | const Vector<Double>& cv = _fitter->compuCovariance().diagonal(); |
1827 1789 | pars.resize(5); |
1828 1790 | pars(0) = solution(1); |
1829 1791 | pars(1) = sqrt(cv(1)); |
1830 1792 | pars(2) = solution(0); |
1831 1793 | pars(3) = sqrt(cv(0)); |
1832 - | pars(4) = itsFitterPtr->chiSquare(); |
1794 + | pars(4) = _fitter->chiSquare(); |
1833 1795 | // |
1834 1796 | return true; |
1835 1797 | } |
1836 1798 | |
1837 1799 | |
1838 - | String ImagePolarimetry::stokesName (ImagePolarimetry::StokesTypes index) const |
1800 + | String ImagePolarimetry::_stokesName (ImagePolarimetry::StokesTypes index) const |
1839 1801 | { |
1840 1802 | if (index==ImagePolarimetry::I) { |
1841 1803 | return String("I"); |
1842 1804 | } else if (index==ImagePolarimetry::Q) { |
1843 1805 | return String("Q"); |
1844 1806 | } else if (index==ImagePolarimetry::U) { |
1845 1807 | return String("U"); |
1846 1808 | } else if (index==ImagePolarimetry::V) { |
1847 1809 | return String("V"); |
1848 1810 | } else { |
1849 1811 | return String("??"); |
1850 1812 | } |
1851 1813 | } |
1852 1814 | |
1853 1815 | |
1854 - | Stokes::StokesTypes ImagePolarimetry::stokesType (ImagePolarimetry::StokesTypes index) const |
1816 + | Stokes::StokesTypes ImagePolarimetry::_stokesType (ImagePolarimetry::StokesTypes index) const |
1855 1817 | { |
1856 1818 | if (index==ImagePolarimetry::I) { |
1857 1819 | return Stokes::I; |
1858 1820 | } else if (index==ImagePolarimetry::Q) { |
1859 1821 | return Stokes::Q; |
1860 1822 | } else if (index==ImagePolarimetry::U) { |
1861 1823 | return Stokes::U; |
1862 1824 | } else if (index==ImagePolarimetry::V) { |
1863 1825 | return Stokes::V; |
1864 1826 | } else { |
1865 1827 | return Stokes::Undefined; |
1866 1828 | } |
1867 1829 | } |
1868 1830 | |
1869 1831 | |
1870 - | Float ImagePolarimetry::sigma (ImagePolarimetry::StokesTypes index, Float clip) |
1832 + | Float ImagePolarimetry::_sigma(ImagePolarimetry::StokesTypes index, Float clip) |
1871 1833 | { |
1872 1834 | Float clip2 = abs(clip); |
1873 1835 | if (clip2==0.0) clip2 = 10.0; |
1874 1836 | // |
1875 - | if (clip2 != itsOldClip && itsStokesStatsPtr[index]!=0) { |
1876 - | delete itsStokesStatsPtr[index]; |
1877 - | itsStokesStatsPtr[index] = 0; |
1837 + | if (clip2 != _oldClip && _stokesStats[index]!=0) { |
1838 + | delete _stokesStats[index]; |
1839 + | _stokesStats[index] = 0; |
1878 1840 | } |
1879 - | if (itsStokesStatsPtr[index]==0) { |
1841 + | if (_stokesStats[index]==0) { |
1880 1842 | |
1881 1843 | // Find sigma for all points inside +/- clip-sigma of the mean |
1882 1844 | // More joys of LEL |
1883 1845 | |
1884 - | const ImageInterface<Float>* p = itsStokesPtr[index]; |
1846 + | const ImageInterface<Float>* p = _stokes[index]; |
1885 1847 | LatticeExprNode n1 (*p); |
1886 1848 | LatticeExprNode n2 (n1[abs(n1-mean(n1)) < clip2*stddev(n1)]); |
1887 1849 | LatticeExpr<Float> le(n2); |
1888 1850 | // |
1889 - | itsStokesStatsPtr[index] = new LatticeStatistics<Float>(le, false, false); |
1851 + | _stokesStats[index] = new LatticeStatistics<Float>(le, false, false); |
1890 1852 | } |
1891 1853 | // |
1892 1854 | Array<Float> sigmaA; |
1893 - | itsStokesStatsPtr[index]->getConvertedStatistic(sigmaA, LatticeStatsBase::SIGMA); |
1855 + | _stokesStats[index]->getConvertedStatistic(sigmaA, LatticeStatsBase::SIGMA); |
1894 1856 | if (sigmaA.nelements()==0) { |
1895 1857 | LogIO os(LogOrigin("ImagePolarimetry", "sigma(...)", WHERE)); |
1896 1858 | os << "No good points in clipped determination of the noise " |
1897 - | << "for the Stokes " << stokesName(index) << " image" << LogIO::EXCEPTION; |
1859 + | << "for the Stokes " << _stokesName(index) << " image" << LogIO::EXCEPTION; |
1898 1860 | } |
1899 1861 | // |
1900 - | itsOldClip = clip2; |
1862 + | _oldClip = clip2; |
1901 1863 | return sigmaA(IPosition(1,0)); |
1902 1864 | } |
1903 1865 | |
1904 1866 | |
1905 - | void ImagePolarimetry::subtractProfileMean (ImageInterface<Float>& im, uInt axis) const |
1867 + | void ImagePolarimetry::_subtractProfileMean (ImageInterface<Float>& im, uInt axis) const |
1906 1868 | { |
1907 1869 | const IPosition tileShape = im.niceCursorShape(); |
1908 1870 | TiledLineStepper ts(im.shape(), tileShape, axis); |
1909 1871 | LatticeIterator<Float> it(im, ts); |
1910 1872 | // |
1911 1873 | Float dMean; |
1912 1874 | if (im.isMasked()) { |
1913 1875 | const Lattice<Bool>& mask = im.pixelMask(); |
1914 1876 | for (it.reset(); !it.atEnd(); it++) { |
1915 1877 | const Array<Float>& p = it.cursor(); |
1922 1884 | |
1923 1885 | } else { |
1924 1886 | for (it.reset(); !it.atEnd(); it++) { |
1925 1887 | dMean = mean(it.vectorCursor()); |
1926 1888 | it.rwVectorCursor() -= dMean; |
1927 1889 | } |
1928 1890 | } |
1929 1891 | } |
1930 1892 | |
1931 1893 | |
1932 - | Bool ImagePolarimetry::dealWithMask (Lattice<Bool>*& pMask, ImageInterface<Float>*& pIm, |
1894 + | Bool ImagePolarimetry::_dealWithMask (Lattice<Bool>*& pMask, ImageInterface<Float>*& pIm, |
1933 1895 | LogIO& os, const String& type) const |
1934 1896 | { |
1935 1897 | Bool isMasked = false; |
1936 1898 | if (!pIm->isMasked()) { |
1937 1899 | if (pIm->canDefineRegion()) { |
1938 1900 | pIm->makeMask("mask0", true, true, true, true); |
1939 1901 | isMasked = true; |
1940 1902 | } else { |
1941 1903 | os << LogIO::WARN << "Could not create a mask for the output " << type << " image" << LogIO::POST; |
1942 1904 | } |
1949 1911 | if (!pMask->isWritable()) { |
1950 1912 | isMasked = false; |
1951 1913 | os << LogIO::WARN << "The output " << type << " image has a mask but it's not writable" << LogIO::POST; |
1952 1914 | } |
1953 1915 | } |
1954 1916 | return isMasked; |
1955 1917 | } |
1956 1918 | |
1957 1919 | void ImagePolarimetry::_createBeamsEqMat() { |
1958 1920 | _beamsEqMat.assign(Matrix<Bool>(4, 4, false)); |
1959 - | Bool hasMultiBeams = itsInImagePtr->imageInfo().hasMultipleBeams(); |
1921 + | Bool hasMultiBeams = _image->imageInfo().hasMultipleBeams(); |
1960 1922 | for (uInt i=0; i<4; i++) { |
1961 1923 | for (uInt j=i; j<4; j++) { |
1962 - | if (itsStokesPtr[i] == 0 || itsStokesPtr[j] == 0) { |
1924 + | if (_stokes[i] == 0 || _stokes[j] == 0) { |
1963 1925 | _beamsEqMat(i, j) = false; |
1964 1926 | } |
1965 1927 | else if (i == j) { |
1966 1928 | _beamsEqMat(i, j) = true; |
1967 1929 | } |
1968 1930 | else if (hasMultiBeams) { |
1969 1931 | _beamsEqMat(i, j) = ( |
1970 - | itsStokesPtr[i]->imageInfo().getBeamSet() |
1971 - | == itsStokesPtr[j]->imageInfo().getBeamSet() |
1932 + | _stokes[i]->imageInfo().getBeamSet() |
1933 + | == _stokes[j]->imageInfo().getBeamSet() |
1972 1934 | ); |
1973 1935 | } |
1974 1936 | else { |
1975 1937 | _beamsEqMat(i, j) = true; |
1976 1938 | } |
1977 1939 | _beamsEqMat(j, i) = _beamsEqMat(i, j); |
1978 1940 | |
1979 1941 | } |
1980 1942 | } |
1981 1943 | } |
2002 1964 | ); |
2003 1965 | } |
2004 1966 | else { |
2005 1967 | return false; |
2006 1968 | } |
2007 1969 | } |
2008 1970 | } |
2009 1971 | } |
2010 1972 | if ( |
2011 1973 | requireChannelEquality |
2012 - | && itsStokesPtr[stokes[0]]->coordinates().hasSpectralAxis() |
2013 - | && itsStokesPtr[stokes[0]]->imageInfo().hasMultipleBeams() |
1974 + | && _stokes[stokes[0]]->coordinates().hasSpectralAxis() |
1975 + | && _stokes[stokes[0]]->imageInfo().hasMultipleBeams() |
2014 1976 | ) { |
2015 - | const Array<GaussianBeam>& beamSet = itsStokesPtr[stokes[0]]->imageInfo().getBeamSet().getBeams(); |
1977 + | const Array<GaussianBeam>& beamSet = _stokes[stokes[0]]->imageInfo().getBeamSet().getBeams(); |
2016 1978 | Array<GaussianBeam>::const_iterator start = beamSet.begin(); |
2017 1979 | start++; |
2018 1980 | for ( |
2019 1981 | Array<GaussianBeam>::const_iterator iter=start; |
2020 1982 | iter!=beamSet.end(); iter++ |
2021 1983 | ) { |
2022 1984 | if (*iter != *(beamSet.begin())) { |
2023 1985 | if (throws) { |
2024 1986 | throw AipsError( |
2025 1987 | "At least one beam in this image is not equal to all the others along its spectral axis so this computation cannot be performed" |