Line data Source code
1 : //# Copyright (C) 1998,1999,2000,2001,2003
2 : //# Associated Universities, Inc. Washington DC, USA.
3 : //#
4 : //# This program is free software; you can redistribute it and/or modify it
5 : //# under the terms of the GNU General Public License as published by the Free
6 : //# Software Foundation; either version 2 of the License, or (at your option)
7 : //# any later version.
8 : //#
9 : //# This program is distributed in the hope that it will be useful, but WITHOUT
10 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
11 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
12 : //# more details.
13 : //#
14 : //# You should have received a copy of the GNU General Public License along
15 : //# with this program; if not, write to the Free Software Foundation, Inc.,
16 : //# 675 Massachusetts Ave, Cambridge, MA 02139, USA.
17 : //#
18 : //# Correspondence concerning AIPS++ should be addressed as follows:
19 : //# Internet email: aips2-request@nrao.edu.
20 : //# Postal address: AIPS++ Project Office
21 : //# National Radio Astronomy Observatory
22 : //# 520 Edgemont Road
23 : //# Charlottesville, VA 22903-2475 USA
24 : //#
25 : //# $Id: $
26 :
27 : #include <imageanalysis/IO/ImageProfileFitterResults.h>
28 :
29 : #include <casacore/casa/Arrays/ArrayLogical.h>
30 : #include <casacore/casa/Utilities/Precision.h>
31 : #include <casacore/coordinates/Coordinates/DirectionCoordinate.h>
32 : #include <casacore/coordinates/Coordinates/LinearCoordinate.h>
33 : #include <casacore/coordinates/Coordinates/SpectralCoordinate.h>
34 : #include <casacore/coordinates/Coordinates/StokesCoordinate.h>
35 : #include <casacore/images/Images/PagedImage.h>
36 : #include <casacore/scimath/Mathematics/Combinatorics.h>
37 :
38 : #include <imageanalysis/ImageAnalysis/ImageCollapser.h>
39 : #include <imageanalysis/ImageAnalysis/ProfileFitResults.h>
40 : #include <imageanalysis/IO/LogFile.h>
41 :
42 : using namespace casacore;
43 : namespace casa {
44 :
45 : const String ImageProfileFitterResults::_class = "ImageProfileFitterResults";
46 : const String ImageProfileFitterResults::_CONVERGED = "converged";
47 : const String ImageProfileFitterResults::_SUCCEEDED = "succeeded";
48 : const String ImageProfileFitterResults::_VALID = "valid";
49 :
50 : const uInt ImageProfileFitterResults::_nOthers = 2;
51 : const uInt ImageProfileFitterResults::_gsPlane = 0;
52 : const uInt ImageProfileFitterResults::_lsPlane = 1;
53 :
54 :
55 86 : ImageProfileFitterResults::ImageProfileFitterResults(
56 : const std::shared_ptr<LogIO> log, const CoordinateSystem& csysIm,
57 : const Array<std::shared_ptr<ProfileFitResults> >* const &fitters, const SpectralList& nonPolyEstimates,
58 : const std::shared_ptr<const SubImage<Float> > subImage, Int polyOrder, Int fitAxis,
59 : uInt nGaussSinglets, uInt nGaussMultiplets, uInt nLorentzSinglets,
60 : uInt nPLPCoeffs, uInt nLTPCoeffs,
61 : Bool logResults, Bool multiFit, const std::shared_ptr<LogFile> logfile,
62 : const String& xUnit, const String& summaryHeader
63 86 : ) : _logResults(logResults), _multiFit(multiFit),
64 : _doVelocity(
65 86 : fitAxis == subImage->coordinates().spectralAxisNumber()
66 172 : && Quantity(1, xUnit).isConform("m/s")
67 172 : && subImage->coordinates().spectralCoordinate().restFrequency() > 0
68 : ), _xUnit(xUnit), _summaryHeader(summaryHeader),
69 : _nGaussSinglets(nGaussSinglets), _nGaussMultiplets(nGaussMultiplets),
70 : _nLorentzSinglets(nLorentzSinglets), _nPLPCoeffs(nPLPCoeffs),_nLTPCoeffs(nLTPCoeffs),
71 : _fitters(fitters),
72 : _nonPolyEstimates(nonPolyEstimates), _subImage(subImage), _polyOrder(polyOrder),
73 258 : _fitAxis(fitAxis), _logfile(logfile), _log(log), _csysIm(csysIm) {}
74 :
75 86 : ImageProfileFitterResults::~ImageProfileFitterResults() {}
76 :
77 83 : Record ImageProfileFitterResults::getResults() const {
78 83 : return _results;
79 : }
80 :
81 83 : void ImageProfileFitterResults::createResults() {
82 83 : _setResults();
83 83 : _resultsToLog();
84 83 : }
85 :
86 83 : std::unique_ptr<vector<vector<Array<Double> > > > ImageProfileFitterResults::_createPCFArrays() const {
87 : std::unique_ptr<vector<vector<Array<Double> > > > pcfArrays(
88 : new vector<vector<Array<Double> > >(
89 166 : NGSOLMATRICES, vector<Array<Double> >(_nGaussMultiplets+_nOthers)
90 83 : )
91 249 : );
92 83 : uInt nSubcomps = 0;
93 83 : uInt compCount = 0;
94 83 : Double fNAN = casacore::doubleNaN();
95 :
96 166 : Array<Double> blank;
97 166 : IPosition fShape = _fitters->shape();
98 251 : for (uInt i=0; i<_nGaussMultiplets + _nOthers; i++) {
99 168 : if (i == _gsPlane) {
100 : // gaussian singlets go in position 0
101 83 : nSubcomps = _nGaussSinglets;
102 : }
103 85 : else if (i == _lsPlane) {
104 : // lorentzian singlets go in position 1
105 83 : nSubcomps = _nLorentzSinglets;
106 : }
107 : else {
108 : // gaussian multiplets go in positions 2 to 1 + _nGaussianMultiplets
109 0 : while (
110 2 : _nonPolyEstimates[compCount]->getType() != SpectralElement::GMULTIPLET
111 : ) {
112 0 : compCount++;
113 : }
114 2 : nSubcomps = dynamic_cast<const GaussianMultipletSpectralElement*>(
115 2 : _nonPolyEstimates[compCount]
116 2 : )->getGaussians().size();
117 2 : compCount++;
118 : }
119 336 : IPosition blankSize(1, nSubcomps);
120 168 : blankSize.prepend(fShape);
121 168 : blank.resize(blankSize, False);
122 168 : blank = fNAN;
123 1512 : for (uInt k=0; k<NGSOLMATRICES; k++) {
124 1344 : (*pcfArrays)[k][i] = blank.copy();
125 : }
126 : }
127 166 : return pcfArrays;
128 : }
129 :
130 86 : vector<String> ImageProfileFitterResults::logSummary(
131 : uInt nProfiles, uInt nAttempted, uInt nSucceeded, uInt nConverged, uInt nValid
132 : ) {
133 86 : vector<String> ret;
134 86 : *_log << LogOrigin(_class, __func__);
135 172 : ostringstream oss;
136 86 : oss << "Number of profiles = " << nProfiles;
137 172 : String str = oss.str();
138 86 : *_log << LogIO::NORMAL << str << LogIO::POST;
139 86 : _writeLogfile(str + "\n", True, False);
140 86 : ret.push_back(str);
141 86 : oss.str("");
142 86 : oss << "Number of fits attempted = " << nAttempted;
143 86 : str = oss.str();
144 86 : *_log << LogOrigin(_class, __func__);
145 86 : *_log << LogIO::NORMAL << str << LogIO::POST;
146 86 : _writeLogfile(str + "\n", False, False);
147 86 : ret.push_back(str);
148 86 : oss.str("");
149 86 : oss << "Number succeeded = " << nSucceeded;
150 86 : str = oss.str();
151 86 : *_log << LogOrigin(_class, __func__);
152 86 : *_log << LogIO::NORMAL << str << LogIO::POST;
153 86 : _writeLogfile(str + "\n", False, False);
154 86 : ret.push_back(str);
155 86 : oss.str("");
156 86 : oss << "Number converged = " << nConverged;
157 86 : str = oss.str();
158 86 : *_log << LogOrigin(_class, __func__);
159 86 : *_log << LogIO::NORMAL << str << LogIO::POST;
160 86 : _writeLogfile(str + "\n", False, False);
161 86 : ret.push_back(str);
162 86 : oss.str("");
163 86 : oss << "Number valid = " << nValid << endl;
164 86 : str = oss.str();
165 86 : *_log << LogOrigin(_class, __func__);
166 86 : *_log << LogIO::NORMAL << str << LogIO::POST;
167 86 : _writeLogfile(str + "\n", False, False);
168 86 : ret.push_back(str);
169 172 : return ret;
170 : }
171 :
172 83 : Bool ImageProfileFitterResults::_setAxisTypes() {
173 83 : const CoordinateSystem& csysSub = _subImage->coordinates();
174 :
175 83 : const DirectionCoordinate *dcoord = csysSub.hasDirectionCoordinate()
176 83 : ? &csysSub.directionCoordinate() : 0;
177 166 : String directionType = dcoord ? MDirection::showType(dcoord->directionType()) : "";
178 83 : const SpectralCoordinate *spcoord = csysSub.hasSpectralAxis()
179 83 : ? &csysSub.spectralCoordinate() : 0;
180 83 : const StokesCoordinate *polcoord = csysSub.hasPolarizationCoordinate()
181 83 : ? &csysSub.stokesCoordinate() : 0;
182 83 : Vector<String> axisNames = csysSub.worldAxisNames();
183 397 : for (uInt i=0; i<axisNames.size(); i++) {
184 314 : axisNames[i].upcase();
185 : }
186 83 : Bool hasLat = False;
187 83 : Bool hasLong = False;
188 83 : _axisTypes = vector<axisType>(axisNames.size(), NAXISTYPES);
189 397 : for (uInt i=0; i<axisNames.size(); i++) {
190 314 : if ((Int)i != _fitAxis) {
191 231 : if (
192 : dcoord
193 231 : && (axisNames[i].startsWith("RIG") || axisNames[i].startsWith("LONG"))
194 : ) {
195 83 : _axisTypes[i] = LONGITUDE;
196 83 : hasLat = True;
197 : }
198 148 : else if (
199 : dcoord
200 148 : && (axisNames[i].startsWith("DEC") || axisNames[i].startsWith("LAT"))
201 : ) {
202 83 : _axisTypes[i] = LATITUDE;
203 83 : hasLong = True;
204 : }
205 65 : else if (
206 : spcoord
207 65 : && (axisNames[i].startsWith("FREQ") || axisNames[i].startsWith("VEL"))
208 : ) {
209 0 : _axisTypes[i] = FREQUENCY;
210 : }
211 65 : else if (
212 : polcoord
213 65 : && axisNames[i].startsWith("STO")
214 : ) {
215 65 : _axisTypes[i] = POLARIZATION;
216 : }
217 : }
218 : }
219 166 : return hasLat && hasLong;
220 : }
221 :
222 83 : void ImageProfileFitterResults::_marshalFitResults(
223 : Array<Bool>& attemptedArr, Array<Bool>& successArr,
224 : Array<Bool>& convergedArr, Array<Bool>& validArr,
225 : Array<String>& typeMat, Array<Int>& niterArr,
226 : Array<Int>& nCompArr, std::unique_ptr<vector<vector<Array<Double> > > >& pcfArrays,
227 : vector<Array<Double> >& plpArrays, vector<Array<Double> >& ltpArrays,
228 : Bool returnDirection, Array<String>& directionInfo
229 : ) {
230 166 : IPosition inTileShape = _subImage->niceCursorShape();
231 166 : TiledLineStepper stepper (_subImage->shape(), inTileShape, _fitAxis);
232 166 : RO_MaskedLatticeIterator<Float> inIter(*_subImage, stepper);
233 83 : const CoordinateSystem& csysSub = _subImage->coordinates();
234 83 : const DirectionCoordinate *dcoord = csysSub.hasDirectionCoordinate()
235 83 : ? &csysSub.directionCoordinate() : 0;
236 166 : String directionType = dcoord ? MDirection::showType(dcoord->directionType()) : "";
237 83 : const SpectralCoordinate *spcoord = csysSub.hasSpectralAxis()
238 83 : ? &csysSub.spectralCoordinate() : 0;
239 83 : const StokesCoordinate *polcoord = csysSub.hasPolarizationCoordinate()
240 83 : ? &csysSub.stokesCoordinate() : 0;
241 166 : IPosition pixel;
242 83 : Double increment = fabs(_fitAxisIncrement());
243 279830 : for (
244 559743 : inIter.reset(); ! inIter.atEnd(); ++inIter
245 : ) {
246 279830 : IPosition pixel = inIter.position();
247 279830 : std::shared_ptr<const ProfileFitResults> fitter = (*_fitters)(pixel);
248 279830 : if (! fitter) {
249 263075 : continue;
250 : }
251 16755 : attemptedArr(pixel) = fitter->getList().nelements() > 0;
252 16755 : successArr(pixel) = fitter->succeeded();
253 50218 : convergedArr(pixel) = attemptedArr(pixel) && successArr(pixel)
254 33463 : && fitter->converged();
255 16755 : validArr(pixel) = convergedArr(pixel) && fitter->isValid();
256 16755 : _doWorldCoords(
257 : directionInfo, csysSub, pixel, dcoord, spcoord,
258 : polcoord, returnDirection, directionType
259 : );
260 16755 : if (validArr(pixel)) {
261 13275 : _processSolutions(
262 : typeMat, niterArr, nCompArr, pixel, fitter,
263 : pcfArrays, plpArrays, ltpArrays, increment
264 : );
265 : }
266 : }
267 83 : }
268 :
269 13275 : void ImageProfileFitterResults::_processSolutions(
270 : Array<String>& typeMat, Array<Int>& niterArr,
271 : Array<Int>& nCompArr, const IPosition& pixel,
272 : std::shared_ptr<const ProfileFitResults> fitter,
273 : std::unique_ptr<vector<vector<Array<Double> > > >& pcfArrays,
274 : vector<Array<Double> >& plpArrays, vector<Array<Double> >& ltpArrays,
275 : Double increment
276 : ) {
277 : // mask(pixel) = anyTrue(inIter.getMask());
278 13275 : niterArr(pixel) = (Int)fitter->getNumberIterations();
279 13275 : nCompArr(pixel) = (Int)fitter->getList().nelements();
280 26550 : SpectralList solutions = fitter->getList();
281 13275 : uInt gCount = 0;
282 13275 : uInt gmCount = 0;
283 13275 : uInt lseCount = 0;
284 13275 : uInt nSolutions = solutions.nelements();
285 39525 : for (uInt i=0; i<nSolutions; i++) {
286 26250 : SpectralElement::Types type = solutions[i]->getType();
287 52500 : IPosition tPos(1, i);
288 26250 : tPos.prepend(pixel);
289 26250 : typeMat(tPos) = SpectralElement::fromType(type);
290 26250 : switch (type) {
291 290 : case SpectralElement::POLYNOMIAL:
292 290 : break;
293 25736 : case SpectralElement::GAUSSIAN:
294 : // allow fall through because gaussians and lorentzians use common code
295 : case SpectralElement::LORENTZIAN:
296 : {
297 0 : const PCFSpectralElement *pcf = dynamic_cast<
298 : const PCFSpectralElement*
299 25736 : >(solutions[i]);
300 25736 : uInt plane = _lsPlane;
301 25736 : uInt col = lseCount;
302 : // if block because we allow case fall through
303 25736 : if (type == SpectralElement::LORENTZIAN) {
304 50 : lseCount++;
305 : }
306 : else {
307 25686 : plane = _gsPlane;
308 25686 : col = gCount;
309 25686 : gCount++;
310 : }
311 25736 : _insertPCF(
312 25736 : *pcfArrays, pixel, *pcf, plane, col,
313 : increment
314 : );
315 : }
316 25736 : break;
317 26 : case SpectralElement::GMULTIPLET:
318 : {
319 0 : const GaussianMultipletSpectralElement *gm = dynamic_cast<
320 : const GaussianMultipletSpectralElement*
321 26 : >(solutions[i]);
322 26 : const Vector<GaussianSpectralElement> g(gm->getGaussians());
323 104 : for (uInt k=0; k<g.size(); k++) {
324 78 : _insertPCF(
325 78 : *pcfArrays, pixel, g[k], gmCount+2,
326 : k, increment
327 : );
328 : }
329 26 : gmCount++;
330 : }
331 26 : break;
332 155 : case SpectralElement::POWERLOGPOLY:
333 : {
334 310 : Vector<Double> sols = solutions[i]->get();
335 310 : Vector<Double> errs = solutions[i]->getError();
336 465 : for (uInt j=0; j<_nPLPCoeffs; j++) {
337 : // here
338 310 : IPosition arrIdx(1, j);
339 310 : arrIdx.prepend(pixel);
340 310 : plpArrays[SPXSOL](arrIdx) = sols[j];
341 310 : plpArrays[SPXERR](arrIdx) = errs[j];
342 : }
343 : }
344 155 : break;
345 43 : case SpectralElement::LOGTRANSPOLY:
346 : {
347 86 : auto sols = solutions[i]->get();
348 86 : auto errs = solutions[i]->getError();
349 129 : for (uInt j=0; j<_nLTPCoeffs; j++) {
350 86 : IPosition arrIdx(1, j);
351 86 : arrIdx.prepend(pixel);
352 86 : ltpArrays[SPXSOL](arrIdx) = sols[j];
353 86 : ltpArrays[SPXERR](arrIdx) = errs[j];
354 : }
355 : }
356 43 : break;
357 0 : default:
358 0 : ThrowCc(
359 : "Logic Error: Unhandled Spectral Element type"
360 : );
361 : break;
362 : }
363 : }
364 13275 : }
365 :
366 16755 : void ImageProfileFitterResults::_doWorldCoords(
367 : Array<String>& directionInfo, const CoordinateSystem& csysSub,
368 : const IPosition& pixel, const DirectionCoordinate* const &dcoord,
369 : const SpectralCoordinate * const &spcoord,
370 : const StokesCoordinate* const &polcoord, Bool returnDirection,
371 : const String& directionType
372 : ) {
373 33510 : Vector<Double> world;
374 16755 : if (csysSub.toWorld(world, pixel)) {
375 33510 : String latitude, longitude;
376 83624 : for (uInt i=0; i<world.size(); i++) {
377 : // The Coordinate::format() calls have the potential of modifying the
378 : // input unit so this needs to be reset at the beginning of the loop
379 : // rather than outside the loop
380 133738 : String emptyUnit = "";
381 66869 : if ((Int)i != _fitAxis) {
382 50114 : if (_axisTypes[i] == LONGITUDE) {
383 16755 : longitude = dcoord->format(emptyUnit, Coordinate::DEFAULT, world[i], 0, True, True);
384 33510 : IPosition x(1, LONGITUDE);
385 16755 : x.append(pixel);
386 16755 : _worldCoords(x) = longitude;
387 : }
388 33359 : else if (_axisTypes[i] == LATITUDE) {
389 16755 : latitude = dcoord->format(emptyUnit, Coordinate::DEFAULT, world[i], 1, True, True);
390 33510 : IPosition x(1, LATITUDE);
391 16755 : x.append(pixel);
392 16755 : _worldCoords(x) = latitude;
393 : }
394 16604 : else if (_axisTypes[i] == FREQUENCY) {
395 0 : IPosition x(1, FREQUENCY);
396 0 : x.append(pixel);
397 0 : _worldCoords(x) = spcoord->format(emptyUnit, Coordinate::DEFAULT, world[i], 0, True, True);
398 : }
399 16604 : else if (_axisTypes[i] == POLARIZATION) {
400 16604 : IPosition x(1, POLARIZATION);
401 16604 : x.append(pixel);
402 16604 : _worldCoords(x) = polcoord->format(emptyUnit, Coordinate::DEFAULT, world[i], 0, True, True);
403 : }
404 : }
405 : }
406 16755 : if (returnDirection) {
407 50265 : directionInfo(pixel) = directionType + " "
408 33510 : + longitude + " " + latitude;
409 : }
410 : }
411 : else {
412 0 : ThrowCc(
413 : "Could not convert pixel to world coordinate: "
414 : + csysSub.errorMessage()
415 : );
416 : }
417 16755 : }
418 :
419 14617 : void ImageProfileFitterResults::_writeLogfile(const String& str, Bool open, Bool close) {
420 14617 : if (_logfile.get() != 0) {
421 221 : _logfile->write(str, open, close);
422 : }
423 14617 : }
424 :
425 83 : void ImageProfileFitterResults::_setResults() {
426 166 : LogOrigin logOrigin(_class, __func__);
427 83 : Double fNAN = casacore::doubleNaN();
428 83 : uInt nComps = _nGaussSinglets + _nGaussMultiplets + _nLorentzSinglets;
429 83 : if (_polyOrder >= 0) {
430 5 : nComps++;
431 : }
432 83 : if (_nPLPCoeffs > 0) {
433 11 : nComps++;
434 : }
435 83 : if (_nLTPCoeffs > 0) {
436 5 : nComps++;
437 : }
438 166 : IPosition fitterShape = _fitters->shape();
439 166 : Array<Bool> attemptedArr(fitterShape, False);
440 166 : Array<Bool> successArr(fitterShape, False);
441 166 : Array<Bool> convergedArr(fitterShape, False);
442 166 : Array<Bool> validArr(fitterShape, False);
443 166 : IPosition wcShape(1, (Int)NAXISTYPES);
444 83 : wcShape.append(fitterShape);
445 83 : _worldCoords = Array<String>(wcShape, String(""));
446 166 : Array<Int> niterArr(fitterShape, -1);
447 : // pfcArrays. Zeroth structure (zeroth vector) index corresponds to
448 : // solution type (amp, center, etc). First structure (first vector)
449 : // index corresponds to type of component (Gaussian, Lorentzian, Gaussian
450 : // multiplet number). Second to n-2 structure (first m Array) indices
451 : // correspond location in the _fitters array. Final structure index
452 : // corresponds to the sub component number (eg for multiple singlets or
453 : // for gaussian multiplet components
454 166 : std::unique_ptr<vector<vector<Array<Double> > > > pcfArrays = _createPCFArrays();
455 166 : IPosition bShape(1, max(_nPLPCoeffs, _nLTPCoeffs));
456 83 : bShape.prepend(fitterShape);
457 166 : Array<Double> blank(bShape, fNAN);
458 : // plpArrays and ltpArrays have the solution type as the zeroth (vector)
459 : // index. The next n indices (the first n indices of the Array) correspond to the position in the _fitters
460 : // array. The final index of the structure (also the final index of the Array) corresponds to
461 : // the componenet number being solved for.
462 166 : vector<Array<Double> > plpArrays(NSPXSOLMATRICES);
463 166 : vector<Array<Double> > ltpArrays(NSPXSOLMATRICES);
464 249 : for (uInt i=0; i<NSPXSOLMATRICES; i++) {
465 : // because assignmet of Arrays is by reference, which is horribly confusing
466 166 : if (_nPLPCoeffs > 0) {
467 22 : plpArrays[i] = blank.copy();
468 : }
469 166 : if (_nLTPCoeffs > 0) {
470 10 : ltpArrays[i] = blank.copy();
471 : }
472 : }
473 166 : IPosition typeShape(1, nComps);
474 83 : typeShape.prepend(fitterShape);
475 166 : Array<String> typeArr(typeShape, String("UNDEF"));
476 166 : Array<Int> nCompArr(fitterShape, -1);
477 83 : Bool returnDirection = _setAxisTypes();
478 166 : Array<String> directionInfo(fitterShape);
479 83 : _marshalFitResults(
480 : attemptedArr, successArr,
481 : convergedArr, validArr,
482 : typeArr, niterArr,
483 : nCompArr, pcfArrays, plpArrays, ltpArrays,
484 : returnDirection, directionInfo
485 : );
486 166 : String key;
487 83 : _results.define("attempted", attemptedArr);
488 83 : _results.define(_SUCCEEDED, successArr);
489 83 : _results.define(_CONVERGED, convergedArr);
490 83 : _results.define(_VALID, validArr);
491 83 : _results.define("niter", niterArr);
492 83 : _results.define("ncomps", nCompArr);
493 83 : if (returnDirection) {
494 83 : _results.define("direction", directionInfo);
495 : }
496 83 : _results.define("xUnit", _xUnit);
497 166 : const String yUnit = _subImage->units().getName();
498 83 : _results.define("yUnit", yUnit);
499 83 : _results.define( "type", typeArr);
500 251 : for (uInt i=0; i<_nGaussMultiplets+_nOthers; ++i) {
501 168 : if (i == _gsPlane && _nGaussSinglets == 0) {
502 101 : continue;
503 : }
504 149 : else if (i == _lsPlane && _nLorentzSinglets == 0) {
505 82 : continue;
506 : }
507 134 : Record rec;
508 67 : rec.define("center", (*pcfArrays)[CENTER][i]);
509 67 : rec.define("fwhm", (*pcfArrays)[FWHM][i]);
510 67 : rec.define("amp", (*pcfArrays)[AMP][i]);
511 67 : rec.define("integral", (*pcfArrays)[INTEGRAL][i]);
512 67 : rec.define("centerErr", (*pcfArrays)[CENTERERR][i]);
513 67 : rec.define("fwhmErr", (*pcfArrays)[FWHMERR][i]);
514 67 : rec.define("ampErr", (*pcfArrays)[AMPERR][i]);
515 67 : rec.define("integralErr", (*pcfArrays)[INTEGRALERR][i]);
516 : String description = i == _gsPlane
517 : ? "Gaussian singlets results"
518 : : i == _lsPlane
519 : ? "Lorentzian singlets"
520 136 : : "Gaussian multiplet number " + String::toString(i-1) + " results";
521 67 : rec.define("description", description);
522 67 : _results.defineRecord(_getTag(i), rec);
523 : }
524 83 : if (_nPLPCoeffs > 0) {
525 11 : Record rec;
526 11 : rec.define("solution", plpArrays[SPXSOL]);
527 11 : rec.define("error", plpArrays[SPXERR]);
528 11 : _results.defineRecord("plp", rec);
529 : }
530 83 : if (_nLTPCoeffs > 0) {
531 5 : Record rec;
532 5 : rec.define("solution", ltpArrays[SPXSOL]);
533 5 : rec.define("error", ltpArrays[SPXERR]);
534 5 : _results.defineRecord("ltp", rec);
535 : }
536 83 : }
537 :
538 72 : String ImageProfileFitterResults::_getTag(const uInt i) const {
539 : return i == _gsPlane
540 : ? "gs"
541 : : i == _lsPlane
542 : ? "ls"
543 141 : : "gm" + String::toString(i-2);
544 : }
545 :
546 25814 : void ImageProfileFitterResults::_insertPCF(
547 : vector<vector<Array<Double> > >& pcfArrays,
548 : const IPosition& pixel, const PCFSpectralElement& pcf,
549 : const uInt idx, const uInt col,
550 : const Double increment
551 : ) const {
552 25814 : IPosition x = pixel;
553 25814 : x.append(IPosition(1, col));
554 25814 : pcfArrays[CENTER][idx](x) = _centerWorld(pcf, pixel);
555 25814 : pcfArrays[FWHM][idx](x) = pcf.getFWHM() * increment;
556 25814 : pcfArrays[AMP][idx](x) = pcf.getAmpl();
557 25814 : pcfArrays[CENTERERR][idx](x) = pcf.getCenterErr() * increment;
558 25814 : pcfArrays[FWHMERR][idx](x) = pcf.getFWHMErr() * increment;
559 25814 : pcfArrays[AMPERR][idx](x) = pcf.getAmplErr();
560 25814 : pcfArrays[INTEGRAL][idx](x) = pcf.getIntegral() * increment;
561 25814 : pcfArrays[INTEGRALERR][idx](x) = pcf.getIntegralErr() * increment;
562 25814 : }
563 :
564 0 : Array<Bool> ImageProfileFitterResults::_replicateMask(
565 : const Array<Bool>& array, Int n
566 : ) {
567 0 : uInt axis = array.ndim() - 1;
568 0 : IPosition newShape = array.shape();
569 0 : newShape[axis] = n;
570 0 : Array<Bool> extended(newShape);
571 0 : IPosition begin(newShape.size(), 0);
572 0 : IPosition end = newShape - 1;
573 0 : for (Int j=0; j<n; j++) {
574 0 : begin[axis] = j;
575 0 : end[axis] = j;
576 0 : extended(begin, end) = array;
577 : }
578 0 : return extended;
579 : }
580 :
581 86 : void ImageProfileFitterResults::writeImages(Bool someConverged) const {
582 : Bool writeSolutionImages = (
583 167 : ! _centerName.empty() || ! _centerErrName.empty()
584 81 : || ! _fwhmName.empty() || ! _fwhmErrName.empty()
585 81 : || ! _ampName.empty() || ! _ampErrName.empty()
586 81 : || ! _integralName.empty() || ! _integralErrName.empty()
587 81 : || ! _plpName.empty() || ! _plpErrName.empty()
588 167 : || ! _ltpName.empty() || ! _ltpErrName.empty()
589 86 : );
590 86 : if (
591 86 : ! _multiFit && writeSolutionImages
592 : ) {
593 0 : *_log << LogIO::WARN << "This was not a multi-pixel fit request so solution "
594 0 : << "images will not be written" << LogIO::POST;
595 : }
596 86 : if (
597 86 : _multiFit && writeSolutionImages
598 : ) {
599 10 : if (
600 10 : _nGaussSinglets + _nGaussMultiplets + _nLorentzSinglets + _nPLPCoeffs + _nLTPCoeffs == 0
601 : ) {
602 0 : *_log << LogIO::WARN << "No functions for which solution images are supported were fit "
603 0 : << "so no solution images will be written" << LogIO::POST;
604 : }
605 : else {
606 10 : if (someConverged) {
607 20 : IPosition axes(1, _fitAxis);
608 : ImageCollapser<Float> collapser(
609 10 : _subImage, axes, False, ImageCollapserData::ZERO, String(""), False
610 30 : );
611 : std::shared_ptr<TempImage<Float> > tmp = std::dynamic_pointer_cast<TempImage<Float> >(
612 10 : collapser.collapse()
613 20 : );
614 10 : ThrowIf(! tmp, "Unable to perform dynamic cast");
615 20 : std::shared_ptr<TempImage<Float> > myTemplate(tmp);
616 20 : const String yUnit = _subImage->units().getName();
617 20 : Array<Bool> mask(_fitters->shape(), False);
618 20 : IPosition inTileShape = _subImage->niceCursorShape();
619 20 : TiledLineStepper stepper (_subImage->shape(), inTileShape, _fitAxis);
620 20 : RO_MaskedLatticeIterator<Float> inIter(*_subImage, stepper);
621 277317 : for (
622 277327 : inIter.reset(); ! inIter.atEnd(); ++inIter
623 : ) {
624 277317 : mask(inIter.position()) = anyTrue(inIter.getMask());
625 : }
626 10 : _writeImages(myTemplate->coordinates(), mask, yUnit);
627 : }
628 : else {
629 0 : *_log << LogIO::WARN << "No solutions converged, solution images will not be written"
630 0 : << LogIO::POST;
631 : }
632 : }
633 : }
634 86 : }
635 :
636 10 : void ImageProfileFitterResults::_writeImages(
637 : const CoordinateSystem& xcsys,
638 : const Array<Bool>& mask, const String& yUnit
639 : ) const {
640 20 : map<String, String> mymap;
641 20 : map<String, String> unitmap;
642 10 : mymap["center"] = _centerName;
643 10 : mymap["centerErr"] = _centerErrName;
644 10 : mymap["fwhm"] = _fwhmName;
645 10 : mymap["fwhmErr"] = _fwhmErrName;
646 10 : mymap["amp"] = _ampName;
647 10 : mymap["ampErr"] = _ampErrName;
648 10 : mymap["integral"] = _integralName;
649 10 : mymap["integralErr"] = _integralErrName;
650 10 : unitmap["center"] = _xUnit;
651 10 : unitmap["centerErr"] = _xUnit;
652 10 : unitmap["fwhm"] = _xUnit;
653 10 : unitmap["fwhmErr"] = _xUnit;
654 10 : unitmap["amp"] = yUnit;
655 10 : unitmap["ampErr"] = yUnit;
656 10 : unitmap["integral"] = _xUnit + "." + yUnit;
657 10 : unitmap["integralErr"] = _xUnit + "." + yUnit;
658 31 : for (uInt i=0; i<_nGaussMultiplets+_nOthers; i++) {
659 21 : if (i == _gsPlane && _nGaussSinglets == 0) {
660 16 : continue;
661 : }
662 14 : else if (i == _lsPlane && _nLorentzSinglets == 0) {
663 9 : continue;
664 : }
665 10 : String id = _getTag(i);
666 45 : for (const auto& p : mymap) {
667 80 : String imagename = p.second;
668 : String suffix = i == _gsPlane
669 : ? ""
670 : : i == _lsPlane
671 : ? "_ls"
672 8 : : _nGaussMultiplets <= 1
673 : ? "_gm"
674 120 : : "_gm" + String::toString(i-_nOthers);
675 40 : imagename += suffix;
676 40 : if (! p.second.empty()) {
677 35 : _makeSolutionImages(
678 : imagename, xcsys,
679 70 : _results.asRecord(id).asArrayDouble(p.first),
680 70 : unitmap.find(p.first)->second, mask
681 : );
682 : }
683 : }
684 : }
685 10 : if (
686 3 : (_nPLPCoeffs > 0 && (! _plpName.empty() || ! _plpErrName.empty()))
687 13 : || (_nLTPCoeffs > 0 && (! _ltpName.empty() || ! _ltpErrName.empty()))
688 : ) {
689 10 : String type = _results.isDefined("plp") ? "plp" : "ltp";
690 5 : if (_nPLPCoeffs > 0 && ! _plpName.empty()) {
691 12 : _makeSolutionImages(
692 3 : _plpName, xcsys,
693 6 : _results.asRecord("plp").asArrayDouble("solution"),
694 : "", mask
695 : );
696 : }
697 5 : if (_nPLPCoeffs > 0 && ! _plpErrName.empty()) {
698 12 : _makeSolutionImages(
699 3 : _plpErrName, xcsys,
700 6 : _results.asRecord("plp").asArrayDouble("error"),
701 : "", mask
702 : );
703 : }
704 5 : if (_nLTPCoeffs > 0 && ! _ltpName.empty()) {
705 8 : _makeSolutionImages(
706 2 : _ltpName, xcsys,
707 4 : _results.asRecord("ltp").asArrayDouble("solution"),
708 : "", mask
709 : );
710 : }
711 5 : if (_nLTPCoeffs > 0 && ! _ltpErrName.empty()) {
712 4 : _makeSolutionImages(
713 1 : _ltpErrName, xcsys,
714 2 : _results.asRecord("ltp").asArrayDouble("error"),
715 : "", mask
716 : );
717 : }
718 : }
719 10 : }
720 : /*
721 : Bool ImageProfileFitterResults::_inVelocitySpace() const {
722 : return _fitAxis == _subImage->coordinates().spectralAxisNumber()
723 : && Quantity(1, _xUnit).isConform("m/s")
724 : && _subImage->coordinates().spectralCoordinate().restFrequency() > 0;
725 : }
726 : */
727 :
728 20969 : Double ImageProfileFitterResults::_fitAxisIncrement() const {
729 20969 : if (_doVelocity) {
730 41906 : Vector<Double> pixels(2);
731 20953 : pixels[0] = 0;
732 20953 : pixels[1] = 1;
733 20953 : Vector<Double> velocities(2);
734 20953 : _subImage->coordinates().spectralCoordinate().pixelToVelocity(
735 : velocities, pixels
736 : );
737 20953 : return velocities[1] - velocities[0];
738 : }
739 : else {
740 16 : return _subImage->coordinates().increment()[_fitAxis];
741 : }
742 : }
743 :
744 0 : const Vector<Double> ImageProfileFitterResults::getPixelCenter(uint index) const {
745 0 : Vector<Double> pos;
746 0 : if ( index < _pixelPositions.size()) {
747 0 : pos = _pixelPositions[index];
748 : }
749 0 : return pos;
750 : }
751 :
752 46700 : Double ImageProfileFitterResults::_centerWorld(
753 : const PCFSpectralElement& solution, const IPosition& imPos
754 : ) const {
755 93400 : Vector<Double> pixel(imPos.size());
756 233089 : for (uInt i=0; i<pixel.size(); i++) {
757 186389 : pixel[i] = imPos[i];
758 : }
759 93400 : Vector<Double> world(pixel.size());
760 : // in pixels here
761 46700 : pixel[_fitAxis] = solution.getCenter();
762 46700 : _subImage->coordinates().toWorld(world, pixel);
763 46700 : if (_doVelocity) {
764 : Double velocity;
765 46700 : _subImage->coordinates().spectralCoordinate().frequencyToVelocity(velocity, world(_fitAxis));
766 46700 : return velocity;
767 : }
768 : else {
769 0 : return world(_fitAxis);
770 : }
771 : }
772 :
773 83 : void ImageProfileFitterResults::_resultsToLog() {
774 83 : if (! _logResults && _logfile.get() == 0) {
775 62 : return;
776 : }
777 42 : ostringstream summary;
778 21 : summary << "****** Fit performed at " << Time().toString() << "******" << endl << endl;
779 21 : summary << _summaryHeader;
780 21 : if (_nPLPCoeffs + _nLTPCoeffs == 0) {
781 5 : if (_goodAmpRange.size() == 2) {
782 0 : summary << " --- valid amplitude range: " << _goodAmpRange << endl;
783 : }
784 5 : if (_goodCenterRange.size() == 2) {
785 0 : summary << " --- valid center range: " << _goodCenterRange << endl;
786 : }
787 5 : if (_goodFWHMRange.size() == 2) {
788 0 : summary << " --- valid FWHM range: " << _goodFWHMRange << endl;
789 : }
790 5 : summary << " --- number of Gaussian singlets: " << _nGaussSinglets << endl;
791 5 : summary << " --- number of Gaussian multiplets: " << _nGaussMultiplets << endl;
792 5 : if (_nGaussMultiplets > 0) {
793 2 : for (uInt i=0; i<_nGaussMultiplets; i++) {
794 3 : Array<Double> amp = _results.asRecord("gm" + String::toString(i)).asArrayDouble(AMP);
795 1 : uInt n = amp.shape()[amp.ndim()-1];
796 1 : summary << " --- number of components in Gaussian multiplet "
797 1 : << i << ": " << n << endl;
798 : }
799 : }
800 5 : if (_polyOrder >= 0) {
801 0 : summary << " --- polynomial order: " << _polyOrder << endl;
802 : }
803 : else {
804 5 : summary << " --- no polynomial fit " << endl;
805 : }
806 : }
807 21 : if (_multiFit) {
808 11 : summary << " --- Multiple profiles fit, one per pixel over selected region" << endl;
809 : }
810 : else {
811 10 : summary << " --- One profile fit, averaged over several pixels if necessary" << endl;
812 : }
813 21 : if (_logResults) {
814 18 : *_log << LogIO::NORMAL << summary.str() << LogIO::POST;
815 : }
816 21 : _writeLogfile(summary.str(), False, False);
817 42 : IPosition inTileShape = _subImage->niceCursorShape();
818 42 : TiledLineStepper stepper (_subImage->shape(), inTileShape, _fitAxis);
819 42 : RO_MaskedLatticeIterator<Float> inIter(*_subImage, stepper);
820 42 : CoordinateSystem csysSub = _subImage->coordinates();
821 42 : Vector<Double> worldStart;
822 21 : ThrowIf(
823 : ! csysSub.toWorld(worldStart, inIter.position()),
824 : csysSub.errorMessage()
825 : );
826 42 : Vector<Double> pixStart;
827 21 : ThrowIf(
828 : ! _csysIm.toPixel(pixStart, worldStart),
829 : _csysIm.errorMessage()
830 : );
831 21 : if (_multiFit) {
832 49 : for (uInt i=0; i<pixStart.size(); i++) {
833 38 : pixStart[i] = (Int)std::floor( pixStart[i] + 0.5);
834 : }
835 : }
836 42 : Vector<Double> imPix(pixStart.size());
837 42 : Vector<Double> world;
838 42 : IPosition subimPos;
839 42 : SpectralList solutions;
840 42 : String axisUnit = _csysIm.worldAxisUnits()[_fitAxis];
841 21 : Int pixPrecision = _multiFit ? 0 : 3;
842 21 : _pixelPositions.resize( _fitters->size());
843 21 : uint index = 0;
844 277216 : for (
845 21 : inIter.reset();
846 277237 : ! inIter.atEnd();
847 277216 : inIter++
848 : ) {
849 277216 : subimPos = inIter.position();
850 277216 : const std::shared_ptr<ProfileFitResults> fitter = (*_fitters)(subimPos);
851 277216 : if (! fitter) {
852 263050 : continue;
853 : }
854 14166 : summary.str("");
855 14166 : Int basePrecision = summary.precision(1);
856 14166 : summary.precision(basePrecision);
857 14166 : if (csysSub.toWorld(world, subimPos)) {
858 14166 : summary << "Fit :" << endl;
859 70730 : for (uInt i=0; i<world.size(); i++) {
860 56564 : if ((Int)i != _fitAxis) {
861 84796 : IPosition x(1, _axisTypes[i]);
862 42398 : x.append(subimPos);
863 42398 : if (_axisTypes[i] == LONGITUDE) {
864 14166 : summary << " RA : "
865 14166 : << _worldCoords(x) << endl;
866 : }
867 28232 : else if (_axisTypes[i] == LATITUDE) {
868 14166 : summary << " Dec : "
869 14166 : << _worldCoords(x) << endl;
870 : }
871 14066 : else if (_axisTypes[i] == FREQUENCY) {
872 0 : summary << " Freq : "
873 0 : << _worldCoords(x) << endl;
874 : }
875 14066 : else if (_axisTypes[i] == POLARIZATION) {
876 14066 : summary << " Stokes : " << _worldCoords(x) << endl;
877 : }
878 : }
879 : }
880 : }
881 : else {
882 0 : ThrowCc(csysSub.errorMessage());
883 : }
884 70730 : for (uInt i=0; i<pixStart.size(); i++) {
885 56564 : imPix[i] = pixStart[i] + subimPos[i];
886 : }
887 14166 : _pixelPositions[index] = imPix;
888 14166 : summary.setf(ios::fixed);
889 14166 : summary << setprecision(pixPrecision) << " Pixel : [";
890 70730 : for (uInt i=0; i<imPix.size(); i++) {
891 56564 : if (i == (uInt)_fitAxis) {
892 14166 : summary << " *";
893 : }
894 : else {
895 42398 : summary << imPix[i];
896 : }
897 56564 : if (i != imPix.size()-1) {
898 42398 : summary << ", ";
899 : }
900 : }
901 14166 : summary << "]" << setprecision(basePrecision) << endl;
902 14166 : summary.unsetf(ios::fixed);
903 14166 : Bool attempted = fitter->getList().nelements() > 0;
904 : summary << " Attempted : "
905 14166 : << (attempted ? "YES" : "NO") << endl;
906 14166 : if (attempted) {
907 28238 : String converged = fitter->converged() ? "YES" : "NO";
908 14119 : summary << " Converged : " << converged << endl;
909 14119 : if (fitter->converged()) {
910 10751 : summary << " Iterations : " << fitter->getNumberIterations() << endl;
911 21502 : String valid = fitter->isValid() ? "YES" : "NO";
912 10751 : summary << " Valid : " << valid << endl;
913 10751 : if (fitter->isValid()) {
914 10751 : solutions = fitter->getList();
915 31785 : for (uInt i=0; i<solutions.nelements(); i++) {
916 21034 : SpectralElement::Types type = solutions[i]->getType();
917 21034 : if (solutions.nelements() > 1) {
918 20566 : summary << " Results for component " << i << ":" << endl;
919 : }
920 21034 : switch(type) {
921 20811 : case SpectralElement::GAUSSIAN:
922 : // allow fall through; gaussians and lorentzians use the same
923 : // method for output
924 : case SpectralElement::LORENTZIAN:
925 : {
926 : const PCFSpectralElement *pcf
927 20811 : = dynamic_cast<const PCFSpectralElement*>(solutions[i]);
928 62433 : summary << _pcfToString(
929 62433 : pcf, _csysIm, world.copy(), subimPos
930 20811 : );
931 : }
932 20811 : break;
933 25 : case SpectralElement::GMULTIPLET:
934 : {
935 : const GaussianMultipletSpectralElement *gm
936 25 : = dynamic_cast<const GaussianMultipletSpectralElement*>(solutions[i]);
937 50 : summary << _gaussianMultipletToString(
938 75 : *gm, _csysIm, world.copy(), subimPos
939 25 : );
940 25 : break;
941 : }
942 0 : case SpectralElement::POLYNOMIAL:
943 : {
944 : const PolynomialSpectralElement *p
945 0 : = dynamic_cast<const PolynomialSpectralElement*>(solutions[i]);
946 0 : summary << _polynomialToString(*p, _csysIm, imPix, world);
947 : }
948 0 : break;
949 155 : case SpectralElement::POWERLOGPOLY:
950 : {
951 : const PowerLogPolynomialSpectralElement *p
952 155 : = dynamic_cast<const PowerLogPolynomialSpectralElement*>(solutions[i]);
953 155 : summary << _powerLogPolyToString(*p);
954 : }
955 155 : break;
956 43 : case SpectralElement::LOGTRANSPOLY:
957 : {
958 : const LogTransformedPolynomialSpectralElement *p
959 43 : = dynamic_cast<const LogTransformedPolynomialSpectralElement*>(solutions[i]);
960 43 : summary << _logTransPolyToString(*p);
961 : }
962 43 : break;
963 0 : default:
964 0 : ThrowCc("Logic Error: Unhandled spectral element type");
965 : break;
966 : }
967 : }
968 : }
969 : }
970 : }
971 14166 : if (_logResults) {
972 13977 : *_log << LogIO::NORMAL << summary.str() << endl << LogIO::POST;
973 : }
974 14166 : _writeLogfile(summary.str(), False, False);
975 : }
976 21 : if (_logfile.get() != 0) {
977 5 : _logfile->close();
978 : }
979 : }
980 :
981 125712 : String ImageProfileFitterResults::_elementToString(
982 : const Double value, const Double error,
983 : const String& unit, Bool isFixed
984 : ) const {
985 251424 : Unit myUnit(unit);
986 251424 : Vector<String> unitPrefix;
987 251424 : String outUnit;
988 251424 : Quantity qVal(value, unit);
989 251424 : Quantity qErr(error, unit);
990 125712 : if (myUnit.getValue() == UnitVal::ANGLE) {
991 0 : Vector<String> angUnits(5);
992 0 : angUnits[0] = "deg";
993 0 : angUnits[1] = "arcmin";
994 0 : angUnits[2] = "arcsec";
995 0 : angUnits[3] = "marcsec";
996 0 : angUnits[4] = "uarcsec";
997 0 : for (uInt i=0; i<angUnits.size(); i++) {
998 0 : outUnit = angUnits[i];
999 0 : if (fabs(qVal.getValue(outUnit)) > 1) {
1000 0 : qVal.convert(outUnit);
1001 0 : qErr.convert(outUnit);
1002 0 : break;
1003 : }
1004 : }
1005 : }
1006 : // some optical images have very weird units that start with numbers
1007 125712 : else if (
1008 250949 : unit.empty() || Quantity(1, myUnit).isConform(Quantity(1, "m/s"))
1009 250949 : || isdigit(unit[0])
1010 : ) {
1011 : // do nothing
1012 : }
1013 : else {
1014 166772 : Vector<String> unitPrefix(10);
1015 83386 : unitPrefix[0] = "T";
1016 83386 : unitPrefix[1] = "G";
1017 83386 : unitPrefix[2] = "M";
1018 83386 : unitPrefix[3] = "k";
1019 83386 : unitPrefix[4] = "";
1020 83386 : unitPrefix[5] = "m";
1021 83386 : unitPrefix[6] = "u";
1022 83386 : unitPrefix[7] = "n";
1023 83386 : unitPrefix[8] = "p";
1024 83386 : unitPrefix[9] = "f";
1025 :
1026 539978 : for (auto prefix: unitPrefix) {
1027 456592 : outUnit = prefix + unit;
1028 456592 : if (fabs(qVal.getValue(outUnit)) > 1) {
1029 82526 : qVal.convert(outUnit);
1030 82526 : qErr.convert(outUnit);
1031 82526 : break;
1032 : }
1033 : }
1034 : }
1035 251424 : Vector<Double> valErr(2);
1036 125712 : valErr[0] = qVal.getValue();
1037 125712 : valErr[1] = qErr.getValue();
1038 :
1039 125712 : uInt precision = precisionForValueErrorPairs(valErr, Vector<Double>());
1040 125712 : ostringstream out;
1041 125712 : out << std::fixed << setprecision(precision);
1042 125712 : out << qVal.getValue();
1043 125712 : if (isFixed && qErr.getValue() == 0) {
1044 12 : out << " (fixed)";
1045 : }
1046 : else {
1047 125700 : out << " +/- " << qErr.getValue();
1048 : }
1049 125712 : out << " " << qVal.getUnit();
1050 251424 : return out.str();
1051 : }
1052 :
1053 20886 : String ImageProfileFitterResults::_pcfToString(
1054 : const PCFSpectralElement *const &pcf, const CoordinateSystem& csys,
1055 : const Vector<Double>& world, const IPosition& imPos,
1056 : const Bool showTypeString, const String& indent
1057 : ) const {
1058 41772 : Vector<Double> myWorld = world;
1059 41772 : String yUnit = _subImage->units().getName();
1060 41772 : ostringstream summary;
1061 41772 : Vector<Bool> fixed = pcf->fixed();
1062 20886 : if (showTypeString) {
1063 20811 : summary << indent << " Type : "
1064 20811 : << SpectralElement::fromType(pcf->getType()) << endl;
1065 : }
1066 20886 : summary << indent << " Peak : "
1067 41772 : << _elementToString(
1068 20886 : pcf->getAmpl(), pcf->getAmplErr(), yUnit, fixed[0]
1069 20886 : ) << endl;
1070 20886 : Double center = _centerWorld(
1071 : *pcf, imPos
1072 : );
1073 :
1074 20886 : Double increment = fabs(_fitAxisIncrement());
1075 :
1076 20886 : Double centerErr = pcf->getCenterErr() * increment;
1077 20886 : Double fwhm = pcf->getFWHM() * increment;
1078 20886 : Double fwhmErr = pcf->getFWHMErr() * increment;
1079 :
1080 20886 : Double pCenter = 0;
1081 20886 : Double pCenterErr = 0;
1082 20886 : Double pFWHM = 0;
1083 20886 : Double pFWHMErr = 0;
1084 20886 : Int specCoordIndex = csys.findCoordinate(Coordinate::SPECTRAL);
1085 20886 : Bool convertedCenterToPix = True;
1086 20886 : Bool convertedFWHMToPix = True;
1087 20886 : if (_doVelocity) {
1088 20886 : if (csys.spectralCoordinate(specCoordIndex).velocityToPixel(pCenter, center)) {
1089 : Double nextVel;
1090 20886 : csys.spectralCoordinate(specCoordIndex).pixelToVelocity(nextVel, pCenter+1);
1091 20886 : Double velInc = fabs(center - nextVel);
1092 20886 : pCenterErr = centerErr/velInc;
1093 20886 : pFWHM = abs(fwhm/velInc);
1094 20886 : pFWHMErr = abs(fwhmErr/velInc);
1095 : }
1096 : else {
1097 0 : convertedCenterToPix = False;
1098 0 : convertedFWHMToPix = False;
1099 : }
1100 : }
1101 : else {
1102 0 : Vector<Double> pixel(myWorld.size());
1103 0 : myWorld[_fitAxis] = center;
1104 0 : Double delta = csys.increment()[_fitAxis];
1105 0 : if (csys.toPixel(pixel, myWorld)) {
1106 0 : pCenter = pixel[_fitAxis];
1107 0 : pCenterErr = abs(centerErr/delta);
1108 : }
1109 : else {
1110 0 : convertedCenterToPix = False;
1111 : }
1112 0 : pFWHM = abs(fwhm/delta);
1113 0 : pFWHMErr = abs(fwhmErr/delta);
1114 : }
1115 20886 : summary << indent << " Center : "
1116 20886 : << _elementToString(
1117 20886 : center, centerErr, _xUnit, fixed[1]
1118 20886 : ) << endl;
1119 20886 : if (convertedCenterToPix) {
1120 20886 : summary << indent << " "
1121 41772 : << _elementToString(
1122 20886 : pCenter, pCenterErr, "pixel", fixed[1]
1123 20886 : ) << endl;
1124 : }
1125 : else {
1126 0 : summary << indent << " Could not convert world to pixel for center" << endl;
1127 : }
1128 20886 : summary << indent << " FWHM : "
1129 20886 : << _elementToString(
1130 20886 : fwhm, fwhmErr, _xUnit, fixed[2]
1131 20886 : )
1132 20886 : << endl;
1133 20886 : if (convertedFWHMToPix) {
1134 20886 : summary << indent << " "
1135 41772 : << _elementToString(pFWHM, pFWHMErr, "pixel", fixed[2])
1136 20886 : << endl;
1137 : }
1138 : else {
1139 0 : summary << indent << " Could not convert FWHM to pixel" << endl;
1140 : }
1141 20886 : Double integral = pcf->getIntegral()*increment;
1142 20886 : Double integralErr = pcf->getIntegralErr()*increment;
1143 41772 : String integUnit = (Quantity(1.0 ,yUnit)*Quantity(1.0, _xUnit)).getUnit();
1144 20886 : summary << indent << " Integral : "
1145 41772 : << _elementToString(integral, integralErr, integUnit, fixed[0] && fixed[2])
1146 20886 : << endl;
1147 20886 : if (fwhm/increment <= 3) {
1148 6912 : summary << indent << "WARNING: The FWHM is only " << (fwhm/increment)
1149 : << " times the channel width. Be aware that instrumental channelization effects"
1150 6912 : << " may be important." << endl;
1151 : }
1152 41772 : return summary.str();
1153 : }
1154 :
1155 25 : String ImageProfileFitterResults::_gaussianMultipletToString(
1156 : const GaussianMultipletSpectralElement& gm,
1157 : const CoordinateSystem& csys,
1158 : const Vector<Double>& world, const IPosition& imPos
1159 : ) const {
1160 50 : Vector<GaussianSpectralElement> g(gm.getGaussians());
1161 25 : ostringstream summary;
1162 25 : summary << " Type : GAUSSIAN MULTIPLET" << endl;
1163 100 : for (uInt i=0; i<g.size(); i++) {
1164 75 : summary << " Results for subcomponent "
1165 75 : << i << ":" << endl;
1166 : summary
1167 150 : << _pcfToString(&g[i], csys, world, imPos, False, " ")
1168 75 : << endl;
1169 : }
1170 50 : return summary.str();
1171 : }
1172 :
1173 0 : String ImageProfileFitterResults::_polynomialToString(
1174 : const PolynomialSpectralElement& poly, const CoordinateSystem& csys,
1175 : const Vector<Double>& imPix, const Vector<Double>& world
1176 : ) const {
1177 0 : ostringstream summary;
1178 0 : summary << " Type: POLYNOMIAL" << endl;
1179 0 : Vector<Double> parms, errs;
1180 0 : poly.get(parms);
1181 0 : poly.getError(errs);
1182 0 : uInt n = parms.size();
1183 0 : for (uInt j=0; j<n; j++) {
1184 0 : String unit = _subImage->units().getName();
1185 0 : if (j > 0) {
1186 0 : if (j == 1) {
1187 0 : unit += "/(pixel)";
1188 : }
1189 : else {
1190 0 : unit += "/((pixel)" + String::toString(j) + ")";
1191 : }
1192 : }
1193 0 : summary << " c" << j << " : "
1194 0 : << _elementToString(parms[j], errs[j], unit, false) << endl;
1195 : }
1196 : // coefficients in pixel coordinates
1197 : Double x0;
1198 0 : Double deltaX = _fitAxisIncrement();
1199 0 : if (_doVelocity) {
1200 0 : csys.spectralCoordinate().pixelToVelocity(x0, 0);
1201 : }
1202 : else {
1203 0 : Vector<Double> p0 = imPix;
1204 0 : p0[_fitAxis] = 0;
1205 0 : Vector<Double> world0 = world;
1206 0 : csys.toWorld(world0, p0);
1207 0 : x0 = world0[_fitAxis];
1208 : // deltaX = csys.increment()[_fitAxis];
1209 : }
1210 0 : Vector<Double> pCoeff(n, 0);
1211 0 : Vector<Double> pCoeffErr(n, 0);
1212 0 : for (uInt j=0; j<n; j++) {
1213 0 : Double sumsq = 0;
1214 0 : for (uInt k=j; k<n; k++) {
1215 0 : Double multiplier = Combinatorics::choose(k, k-j)
1216 0 : * casacore::pow(x0, Float(k - j))
1217 0 : * casacore::pow(1/deltaX, Float(k));
1218 0 : if ((k-j) % 2 == 1) {
1219 0 : multiplier *= -1;
1220 : }
1221 0 : pCoeff[j] += multiplier * parms[k];
1222 0 : Double errCoeff = multiplier * errs[k];
1223 0 : sumsq += errCoeff * errCoeff;
1224 : }
1225 0 : pCoeffErr[j] = casacore::sqrt(sumsq);
1226 0 : summary << " c" << j << " : ";
1227 0 : String unit = _subImage->units().getName();
1228 0 : if (j > 0 ) {
1229 0 : if (j == 1) {
1230 0 : unit += "/(" + _xUnit + ")";
1231 : }
1232 : else {
1233 0 : unit += "/((" + _xUnit + ")" + String::toString(j) + ")";
1234 : }
1235 : }
1236 0 : summary << _elementToString(pCoeff[j], pCoeffErr[j], unit, false) << endl;
1237 : }
1238 0 : return summary.str();
1239 : }
1240 :
1241 :
1242 155 : String ImageProfileFitterResults::_powerLogPolyToString(
1243 : const PowerLogPolynomialSpectralElement& plp
1244 : ) const {
1245 310 : ostringstream summary;
1246 155 : summary << " Type : POWER LOGARITHMIC POLYNOMIAL" << endl;
1247 155 : summary << " Function : c0*(x/D)**(c1";
1248 155 : uInt nElements = plp.get().size();
1249 155 : for (uInt i=2; i<nElements; i++) {
1250 0 : summary << " + c" << i << "*ln(x/D)";
1251 0 : if (i > 2) {
1252 0 : summary << "**" << (i - 1);
1253 : }
1254 : }
1255 155 : summary << ")" << endl;
1256 155 : summary << " D : " << _plpDivisor << endl;
1257 310 : Vector<Double> parms = plp.get();
1258 310 : Vector<Double> errs = plp.getError();
1259 155 : Vector<Bool> fixed = plp.fixed();
1260 :
1261 465 : for (uInt j=0; j<parms.size(); j++) {
1262 310 : summary << " c" << j << " : "
1263 310 : << _elementToString(parms[j], errs[j], "", fixed[j]) << endl;
1264 : }
1265 310 : return summary.str();
1266 : }
1267 :
1268 43 : String ImageProfileFitterResults::_logTransPolyToString(
1269 : const LogTransformedPolynomialSpectralElement& ltp
1270 : ) const {
1271 86 : ostringstream summary;
1272 43 : summary << " Type : LOGARITHMIC TRANSFORMED POLYNOMIAL" << endl;
1273 43 : summary << " Function : ln(y) = c0";
1274 43 : uInt nElements = ltp.get().size();
1275 86 : for (uInt i=1; i<nElements; i++) {
1276 43 : summary << " + c" << i << "*ln(x/D)";
1277 43 : if (i > 2) {
1278 0 : summary << "**" << i;
1279 : }
1280 : }
1281 43 : summary << ")" << endl;
1282 43 : summary << " D : " << _plpDivisor << endl;
1283 86 : Vector<Double> parms = ltp.get();
1284 86 : Vector<Double> errs = ltp.getError();
1285 43 : Vector<Bool> fixed = ltp.fixed();
1286 :
1287 129 : for (uInt j=0; j<parms.size(); j++) {
1288 86 : summary << " c" << j << " : "
1289 86 : << _elementToString(parms[j], errs[j], "", fixed[j]) << endl;
1290 : }
1291 86 : return summary.str();
1292 : }
1293 :
1294 44 : void ImageProfileFitterResults::_makeSolutionImages(
1295 : const String& name, const CoordinateSystem& csys,
1296 : const Array<Double>& values, const String& unit,
1297 : const Array<Bool>& mask
1298 : ) {
1299 88 : auto valuesShape = values.shape();
1300 88 : Vector<Float> dataCopy(values.size());
1301 88 : Vector<Double>::const_iterator iter;
1302 : // isNaN(Array<Double>&) works, isNaN(Array<Float>&) gives spurious results
1303 132 : Array<Bool> nanInfMask = ! (isNaN(values) || isInf(values));
1304 88 : Vector<Float>::iterator jiter = dataCopy.begin();
1305 1139760 : for (iter=values.begin(); iter!=values.end(); ++iter, ++jiter) {
1306 1139716 : *jiter = (Float)*iter;
1307 : }
1308 88 : auto dc = dataCopy.reform(valuesShape);
1309 44 : uInt nImages = valuesShape.last();
1310 88 : auto imageShape = valuesShape.getFirst(valuesShape.size() - 1);
1311 88 : IPosition start(values.ndim(), 0);
1312 88 : IPosition end = valuesShape - 1;
1313 44 : end[end.size() - 1] = 0;
1314 140 : for (uInt i=0; i<nImages; ++i) {
1315 : auto myname = nImages == 0
1316 : ? name
1317 288 : : name + "_" + String::toString(i);
1318 192 : PagedImage<Float> image(imageShape, csys, myname);
1319 96 : start[start.size() - 1] = i;
1320 96 : end[end.size() - 1] = i;
1321 192 : auto imageVals = dc(start, end);
1322 192 : Array<Double> doubleValues = values(start, end);
1323 : // isNaN(Array<Double>&) works, isNaN(Array<Float>&) gives spurious results
1324 : // Its important to use isInf on the float values not the double values
1325 192 : Array<Bool> nanInfMask = ! (isNaN(doubleValues) || isInf(imageVals));
1326 : // remove the last axis
1327 96 : imageVals.removeDegenerate(values.ndim() -1);
1328 96 : image.put(imageVals);
1329 96 : nanInfMask.removeDegenerate(values.ndim() -1);
1330 96 : Bool hasPixMask = ! allTrue(mask);
1331 96 : Bool hasNanMask = ! allTrue(nanInfMask);
1332 96 : if (hasNanMask || hasPixMask) {
1333 68 : Array<Bool> resMask(imageShape);
1334 : String maskName = image.makeUniqueRegionName(
1335 34 : String("mask"), 0
1336 68 : );
1337 34 : image.makeMask(maskName, True, True, False);
1338 34 : if (hasPixMask) {
1339 20 : resMask = mask;
1340 20 : if (hasNanMask) {
1341 20 : resMask = resMask && nanInfMask;
1342 : }
1343 : }
1344 : else {
1345 14 : resMask = nanInfMask;
1346 : }
1347 34 : (&image.pixelMask())->put(resMask);
1348 : }
1349 96 : image.setUnits(Unit(unit));
1350 : }
1351 44 : }
1352 :
1353 : }
|
