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 0 : 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 0 : ) : _logResults(logResults), _multiFit(multiFit),
64 : _doVelocity(
65 0 : fitAxis == subImage->coordinates().spectralAxisNumber()
66 0 : && Quantity(1, xUnit).isConform("m/s")
67 0 : && 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 0 : _fitAxis(fitAxis), _logfile(logfile), _log(log), _csysIm(csysIm) {}
74 :
75 0 : ImageProfileFitterResults::~ImageProfileFitterResults() {}
76 :
77 0 : Record ImageProfileFitterResults::getResults() const {
78 0 : return _results;
79 : }
80 :
81 0 : void ImageProfileFitterResults::createResults() {
82 0 : _setResults();
83 0 : _resultsToLog();
84 0 : }
85 :
86 0 : 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 0 : NGSOLMATRICES, vector<Array<Double> >(_nGaussMultiplets+_nOthers)
90 0 : )
91 0 : );
92 0 : uInt nSubcomps = 0;
93 0 : uInt compCount = 0;
94 0 : Double fNAN = casacore::doubleNaN();
95 :
96 0 : Array<Double> blank;
97 0 : IPosition fShape = _fitters->shape();
98 0 : for (uInt i=0; i<_nGaussMultiplets + _nOthers; i++) {
99 0 : if (i == _gsPlane) {
100 : // gaussian singlets go in position 0
101 0 : nSubcomps = _nGaussSinglets;
102 : }
103 0 : else if (i == _lsPlane) {
104 : // lorentzian singlets go in position 1
105 0 : nSubcomps = _nLorentzSinglets;
106 : }
107 : else {
108 : // gaussian multiplets go in positions 2 to 1 + _nGaussianMultiplets
109 0 : while (
110 0 : _nonPolyEstimates[compCount]->getType() != SpectralElement::GMULTIPLET
111 : ) {
112 0 : compCount++;
113 : }
114 0 : nSubcomps = dynamic_cast<const GaussianMultipletSpectralElement*>(
115 0 : _nonPolyEstimates[compCount]
116 0 : )->getGaussians().size();
117 0 : compCount++;
118 : }
119 0 : IPosition blankSize(1, nSubcomps);
120 0 : blankSize.prepend(fShape);
121 0 : blank.resize(blankSize, False);
122 0 : blank = fNAN;
123 0 : for (uInt k=0; k<NGSOLMATRICES; k++) {
124 0 : (*pcfArrays)[k][i] = blank.copy();
125 : }
126 : }
127 0 : return pcfArrays;
128 : }
129 :
130 0 : vector<String> ImageProfileFitterResults::logSummary(
131 : uInt nProfiles, uInt nAttempted, uInt nSucceeded, uInt nConverged, uInt nValid
132 : ) {
133 0 : vector<String> ret;
134 0 : *_log << LogOrigin(_class, __func__);
135 0 : ostringstream oss;
136 0 : oss << "Number of profiles = " << nProfiles;
137 0 : String str = oss.str();
138 0 : *_log << LogIO::NORMAL << str << LogIO::POST;
139 0 : _writeLogfile(str + "\n", True, False);
140 0 : ret.push_back(str);
141 0 : oss.str("");
142 0 : oss << "Number of fits attempted = " << nAttempted;
143 0 : str = oss.str();
144 0 : *_log << LogOrigin(_class, __func__);
145 0 : *_log << LogIO::NORMAL << str << LogIO::POST;
146 0 : _writeLogfile(str + "\n", False, False);
147 0 : ret.push_back(str);
148 0 : oss.str("");
149 0 : oss << "Number succeeded = " << nSucceeded;
150 0 : str = oss.str();
151 0 : *_log << LogOrigin(_class, __func__);
152 0 : *_log << LogIO::NORMAL << str << LogIO::POST;
153 0 : _writeLogfile(str + "\n", False, False);
154 0 : ret.push_back(str);
155 0 : oss.str("");
156 0 : oss << "Number converged = " << nConverged;
157 0 : str = oss.str();
158 0 : *_log << LogOrigin(_class, __func__);
159 0 : *_log << LogIO::NORMAL << str << LogIO::POST;
160 0 : _writeLogfile(str + "\n", False, False);
161 0 : ret.push_back(str);
162 0 : oss.str("");
163 0 : oss << "Number valid = " << nValid << endl;
164 0 : str = oss.str();
165 0 : *_log << LogOrigin(_class, __func__);
166 0 : *_log << LogIO::NORMAL << str << LogIO::POST;
167 0 : _writeLogfile(str + "\n", False, False);
168 0 : ret.push_back(str);
169 0 : return ret;
170 : }
171 :
172 0 : Bool ImageProfileFitterResults::_setAxisTypes() {
173 0 : const CoordinateSystem& csysSub = _subImage->coordinates();
174 :
175 0 : const DirectionCoordinate *dcoord = csysSub.hasDirectionCoordinate()
176 0 : ? &csysSub.directionCoordinate() : 0;
177 0 : String directionType = dcoord ? MDirection::showType(dcoord->directionType()) : "";
178 0 : const SpectralCoordinate *spcoord = csysSub.hasSpectralAxis()
179 0 : ? &csysSub.spectralCoordinate() : 0;
180 0 : const StokesCoordinate *polcoord = csysSub.hasPolarizationCoordinate()
181 0 : ? &csysSub.stokesCoordinate() : 0;
182 0 : Vector<String> axisNames = csysSub.worldAxisNames();
183 0 : for (uInt i=0; i<axisNames.size(); i++) {
184 0 : axisNames[i].upcase();
185 : }
186 0 : Bool hasLat = False;
187 0 : Bool hasLong = False;
188 0 : _axisTypes = vector<axisType>(axisNames.size(), NAXISTYPES);
189 0 : for (uInt i=0; i<axisNames.size(); i++) {
190 0 : if ((Int)i != _fitAxis) {
191 0 : if (
192 : dcoord
193 0 : && (axisNames[i].startsWith("RIG") || axisNames[i].startsWith("LONG"))
194 : ) {
195 0 : _axisTypes[i] = LONGITUDE;
196 0 : hasLat = True;
197 : }
198 0 : else if (
199 : dcoord
200 0 : && (axisNames[i].startsWith("DEC") || axisNames[i].startsWith("LAT"))
201 : ) {
202 0 : _axisTypes[i] = LATITUDE;
203 0 : hasLong = True;
204 : }
205 0 : else if (
206 : spcoord
207 0 : && (axisNames[i].startsWith("FREQ") || axisNames[i].startsWith("VEL"))
208 : ) {
209 0 : _axisTypes[i] = FREQUENCY;
210 : }
211 0 : else if (
212 : polcoord
213 0 : && axisNames[i].startsWith("STO")
214 : ) {
215 0 : _axisTypes[i] = POLARIZATION;
216 : }
217 : }
218 : }
219 0 : return hasLat && hasLong;
220 : }
221 :
222 0 : 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 0 : IPosition inTileShape = _subImage->niceCursorShape();
231 0 : TiledLineStepper stepper (_subImage->shape(), inTileShape, _fitAxis);
232 0 : RO_MaskedLatticeIterator<Float> inIter(*_subImage, stepper);
233 0 : const CoordinateSystem& csysSub = _subImage->coordinates();
234 0 : const DirectionCoordinate *dcoord = csysSub.hasDirectionCoordinate()
235 0 : ? &csysSub.directionCoordinate() : 0;
236 0 : String directionType = dcoord ? MDirection::showType(dcoord->directionType()) : "";
237 0 : const SpectralCoordinate *spcoord = csysSub.hasSpectralAxis()
238 0 : ? &csysSub.spectralCoordinate() : 0;
239 0 : const StokesCoordinate *polcoord = csysSub.hasPolarizationCoordinate()
240 0 : ? &csysSub.stokesCoordinate() : 0;
241 0 : IPosition pixel;
242 0 : Double increment = fabs(_fitAxisIncrement());
243 0 : for (
244 0 : inIter.reset(); ! inIter.atEnd(); ++inIter
245 : ) {
246 0 : IPosition pixel = inIter.position();
247 0 : std::shared_ptr<const ProfileFitResults> fitter = (*_fitters)(pixel);
248 0 : if (! fitter) {
249 0 : continue;
250 : }
251 0 : attemptedArr(pixel) = fitter->getList().nelements() > 0;
252 0 : successArr(pixel) = fitter->succeeded();
253 0 : convergedArr(pixel) = attemptedArr(pixel) && successArr(pixel)
254 0 : && fitter->converged();
255 0 : validArr(pixel) = convergedArr(pixel) && fitter->isValid();
256 0 : _doWorldCoords(
257 : directionInfo, csysSub, pixel, dcoord, spcoord,
258 : polcoord, returnDirection, directionType
259 : );
260 0 : if (validArr(pixel)) {
261 0 : _processSolutions(
262 : typeMat, niterArr, nCompArr, pixel, fitter,
263 : pcfArrays, plpArrays, ltpArrays, increment
264 : );
265 : }
266 : }
267 0 : }
268 :
269 0 : 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 0 : niterArr(pixel) = (Int)fitter->getNumberIterations();
279 0 : nCompArr(pixel) = (Int)fitter->getList().nelements();
280 0 : SpectralList solutions = fitter->getList();
281 0 : uInt gCount = 0;
282 0 : uInt gmCount = 0;
283 0 : uInt lseCount = 0;
284 0 : uInt nSolutions = solutions.nelements();
285 0 : for (uInt i=0; i<nSolutions; i++) {
286 0 : SpectralElement::Types type = solutions[i]->getType();
287 0 : IPosition tPos(1, i);
288 0 : tPos.prepend(pixel);
289 0 : typeMat(tPos) = SpectralElement::fromType(type);
290 0 : switch (type) {
291 0 : case SpectralElement::POLYNOMIAL:
292 0 : break;
293 0 : 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 0 : >(solutions[i]);
300 0 : uInt plane = _lsPlane;
301 0 : uInt col = lseCount;
302 : // if block because we allow case fall through
303 0 : if (type == SpectralElement::LORENTZIAN) {
304 0 : lseCount++;
305 : }
306 : else {
307 0 : plane = _gsPlane;
308 0 : col = gCount;
309 0 : gCount++;
310 : }
311 0 : _insertPCF(
312 0 : *pcfArrays, pixel, *pcf, plane, col,
313 : increment
314 : );
315 : }
316 0 : break;
317 0 : case SpectralElement::GMULTIPLET:
318 : {
319 0 : const GaussianMultipletSpectralElement *gm = dynamic_cast<
320 : const GaussianMultipletSpectralElement*
321 0 : >(solutions[i]);
322 0 : const Vector<GaussianSpectralElement> g(gm->getGaussians());
323 0 : for (uInt k=0; k<g.size(); k++) {
324 0 : _insertPCF(
325 0 : *pcfArrays, pixel, g[k], gmCount+2,
326 : k, increment
327 : );
328 : }
329 0 : gmCount++;
330 : }
331 0 : break;
332 0 : case SpectralElement::POWERLOGPOLY:
333 : {
334 0 : Vector<Double> sols = solutions[i]->get();
335 0 : Vector<Double> errs = solutions[i]->getError();
336 0 : for (uInt j=0; j<_nPLPCoeffs; j++) {
337 : // here
338 0 : IPosition arrIdx(1, j);
339 0 : arrIdx.prepend(pixel);
340 0 : plpArrays[SPXSOL](arrIdx) = sols[j];
341 0 : plpArrays[SPXERR](arrIdx) = errs[j];
342 : }
343 : }
344 0 : break;
345 0 : case SpectralElement::LOGTRANSPOLY:
346 : {
347 0 : auto sols = solutions[i]->get();
348 0 : auto errs = solutions[i]->getError();
349 0 : for (uInt j=0; j<_nLTPCoeffs; j++) {
350 0 : IPosition arrIdx(1, j);
351 0 : arrIdx.prepend(pixel);
352 0 : ltpArrays[SPXSOL](arrIdx) = sols[j];
353 0 : ltpArrays[SPXERR](arrIdx) = errs[j];
354 : }
355 : }
356 0 : break;
357 0 : default:
358 0 : ThrowCc(
359 : "Logic Error: Unhandled Spectral Element type"
360 : );
361 : break;
362 : }
363 : }
364 0 : }
365 :
366 0 : 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 0 : Vector<Double> world;
374 0 : if (csysSub.toWorld(world, pixel)) {
375 0 : String latitude, longitude;
376 0 : 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 0 : String emptyUnit = "";
381 0 : if ((Int)i != _fitAxis) {
382 0 : if (_axisTypes[i] == LONGITUDE) {
383 0 : longitude = dcoord->format(emptyUnit, Coordinate::DEFAULT, world[i], 0, True, True);
384 0 : IPosition x(1, LONGITUDE);
385 0 : x.append(pixel);
386 0 : _worldCoords(x) = longitude;
387 : }
388 0 : else if (_axisTypes[i] == LATITUDE) {
389 0 : latitude = dcoord->format(emptyUnit, Coordinate::DEFAULT, world[i], 1, True, True);
390 0 : IPosition x(1, LATITUDE);
391 0 : x.append(pixel);
392 0 : _worldCoords(x) = latitude;
393 : }
394 0 : 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 0 : else if (_axisTypes[i] == POLARIZATION) {
400 0 : IPosition x(1, POLARIZATION);
401 0 : x.append(pixel);
402 0 : _worldCoords(x) = polcoord->format(emptyUnit, Coordinate::DEFAULT, world[i], 0, True, True);
403 : }
404 : }
405 : }
406 0 : if (returnDirection) {
407 0 : directionInfo(pixel) = directionType + " "
408 0 : + longitude + " " + latitude;
409 : }
410 : }
411 : else {
412 0 : ThrowCc(
413 : "Could not convert pixel to world coordinate: "
414 : + csysSub.errorMessage()
415 : );
416 : }
417 0 : }
418 :
419 0 : void ImageProfileFitterResults::_writeLogfile(const String& str, Bool open, Bool close) {
420 0 : if (_logfile.get() != 0) {
421 0 : _logfile->write(str, open, close);
422 : }
423 0 : }
424 :
425 0 : void ImageProfileFitterResults::_setResults() {
426 0 : LogOrigin logOrigin(_class, __func__);
427 0 : Double fNAN = casacore::doubleNaN();
428 0 : uInt nComps = _nGaussSinglets + _nGaussMultiplets + _nLorentzSinglets;
429 0 : if (_polyOrder >= 0) {
430 0 : nComps++;
431 : }
432 0 : if (_nPLPCoeffs > 0) {
433 0 : nComps++;
434 : }
435 0 : if (_nLTPCoeffs > 0) {
436 0 : nComps++;
437 : }
438 0 : IPosition fitterShape = _fitters->shape();
439 0 : Array<Bool> attemptedArr(fitterShape, False);
440 0 : Array<Bool> successArr(fitterShape, False);
441 0 : Array<Bool> convergedArr(fitterShape, False);
442 0 : Array<Bool> validArr(fitterShape, False);
443 0 : IPosition wcShape(1, (Int)NAXISTYPES);
444 0 : wcShape.append(fitterShape);
445 0 : _worldCoords = Array<String>(wcShape, String(""));
446 0 : 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 0 : std::unique_ptr<vector<vector<Array<Double> > > > pcfArrays = _createPCFArrays();
455 0 : IPosition bShape(1, max(_nPLPCoeffs, _nLTPCoeffs));
456 0 : bShape.prepend(fitterShape);
457 0 : 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 0 : vector<Array<Double> > plpArrays(NSPXSOLMATRICES);
463 0 : vector<Array<Double> > ltpArrays(NSPXSOLMATRICES);
464 0 : for (uInt i=0; i<NSPXSOLMATRICES; i++) {
465 : // because assignmet of Arrays is by reference, which is horribly confusing
466 0 : if (_nPLPCoeffs > 0) {
467 0 : plpArrays[i] = blank.copy();
468 : }
469 0 : if (_nLTPCoeffs > 0) {
470 0 : ltpArrays[i] = blank.copy();
471 : }
472 : }
473 0 : IPosition typeShape(1, nComps);
474 0 : typeShape.prepend(fitterShape);
475 0 : Array<String> typeArr(typeShape, String("UNDEF"));
476 0 : Array<Int> nCompArr(fitterShape, -1);
477 0 : Bool returnDirection = _setAxisTypes();
478 0 : Array<String> directionInfo(fitterShape);
479 0 : _marshalFitResults(
480 : attemptedArr, successArr,
481 : convergedArr, validArr,
482 : typeArr, niterArr,
483 : nCompArr, pcfArrays, plpArrays, ltpArrays,
484 : returnDirection, directionInfo
485 : );
486 0 : String key;
487 0 : _results.define("attempted", attemptedArr);
488 0 : _results.define(_SUCCEEDED, successArr);
489 0 : _results.define(_CONVERGED, convergedArr);
490 0 : _results.define(_VALID, validArr);
491 0 : _results.define("niter", niterArr);
492 0 : _results.define("ncomps", nCompArr);
493 0 : if (returnDirection) {
494 0 : _results.define("direction", directionInfo);
495 : }
496 0 : _results.define("xUnit", _xUnit);
497 0 : const String yUnit = _subImage->units().getName();
498 0 : _results.define("yUnit", yUnit);
499 0 : _results.define( "type", typeArr);
500 0 : for (uInt i=0; i<_nGaussMultiplets+_nOthers; ++i) {
501 0 : if (i == _gsPlane && _nGaussSinglets == 0) {
502 0 : continue;
503 : }
504 0 : else if (i == _lsPlane && _nLorentzSinglets == 0) {
505 0 : continue;
506 : }
507 0 : Record rec;
508 0 : rec.define("center", (*pcfArrays)[CENTER][i]);
509 0 : rec.define("fwhm", (*pcfArrays)[FWHM][i]);
510 0 : rec.define("amp", (*pcfArrays)[AMP][i]);
511 0 : rec.define("integral", (*pcfArrays)[INTEGRAL][i]);
512 0 : rec.define("centerErr", (*pcfArrays)[CENTERERR][i]);
513 0 : rec.define("fwhmErr", (*pcfArrays)[FWHMERR][i]);
514 0 : rec.define("ampErr", (*pcfArrays)[AMPERR][i]);
515 0 : rec.define("integralErr", (*pcfArrays)[INTEGRALERR][i]);
516 : String description = i == _gsPlane
517 : ? "Gaussian singlets results"
518 : : i == _lsPlane
519 : ? "Lorentzian singlets"
520 0 : : "Gaussian multiplet number " + String::toString(i-1) + " results";
521 0 : rec.define("description", description);
522 0 : _results.defineRecord(_getTag(i), rec);
523 : }
524 0 : if (_nPLPCoeffs > 0) {
525 0 : Record rec;
526 0 : rec.define("solution", plpArrays[SPXSOL]);
527 0 : rec.define("error", plpArrays[SPXERR]);
528 0 : _results.defineRecord("plp", rec);
529 : }
530 0 : if (_nLTPCoeffs > 0) {
531 0 : Record rec;
532 0 : rec.define("solution", ltpArrays[SPXSOL]);
533 0 : rec.define("error", ltpArrays[SPXERR]);
534 0 : _results.defineRecord("ltp", rec);
535 : }
536 0 : }
537 :
538 0 : String ImageProfileFitterResults::_getTag(const uInt i) const {
539 : return i == _gsPlane
540 : ? "gs"
541 : : i == _lsPlane
542 : ? "ls"
543 0 : : "gm" + String::toString(i-2);
544 : }
545 :
546 0 : 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 0 : IPosition x = pixel;
553 0 : x.append(IPosition(1, col));
554 0 : pcfArrays[CENTER][idx](x) = _centerWorld(pcf, pixel);
555 0 : pcfArrays[FWHM][idx](x) = pcf.getFWHM() * increment;
556 0 : pcfArrays[AMP][idx](x) = pcf.getAmpl();
557 0 : pcfArrays[CENTERERR][idx](x) = pcf.getCenterErr() * increment;
558 0 : pcfArrays[FWHMERR][idx](x) = pcf.getFWHMErr() * increment;
559 0 : pcfArrays[AMPERR][idx](x) = pcf.getAmplErr();
560 0 : pcfArrays[INTEGRAL][idx](x) = pcf.getIntegral() * increment;
561 0 : pcfArrays[INTEGRALERR][idx](x) = pcf.getIntegralErr() * increment;
562 0 : }
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 0 : void ImageProfileFitterResults::writeImages(Bool someConverged) const {
582 : Bool writeSolutionImages = (
583 0 : ! _centerName.empty() || ! _centerErrName.empty()
584 0 : || ! _fwhmName.empty() || ! _fwhmErrName.empty()
585 0 : || ! _ampName.empty() || ! _ampErrName.empty()
586 0 : || ! _integralName.empty() || ! _integralErrName.empty()
587 0 : || ! _plpName.empty() || ! _plpErrName.empty()
588 0 : || ! _ltpName.empty() || ! _ltpErrName.empty()
589 0 : );
590 0 : if (
591 0 : ! _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 0 : if (
597 0 : _multiFit && writeSolutionImages
598 : ) {
599 0 : if (
600 0 : _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 0 : if (someConverged) {
607 0 : IPosition axes(1, _fitAxis);
608 : ImageCollapser<Float> collapser(
609 0 : _subImage, axes, False, ImageCollapserData::ZERO, String(""), False
610 0 : );
611 : std::shared_ptr<TempImage<Float> > tmp = std::dynamic_pointer_cast<TempImage<Float> >(
612 0 : collapser.collapse()
613 0 : );
614 0 : ThrowIf(! tmp, "Unable to perform dynamic cast");
615 0 : std::shared_ptr<TempImage<Float> > myTemplate(tmp);
616 0 : const String yUnit = _subImage->units().getName();
617 0 : Array<Bool> mask(_fitters->shape(), False);
618 0 : IPosition inTileShape = _subImage->niceCursorShape();
619 0 : TiledLineStepper stepper (_subImage->shape(), inTileShape, _fitAxis);
620 0 : RO_MaskedLatticeIterator<Float> inIter(*_subImage, stepper);
621 0 : for (
622 0 : inIter.reset(); ! inIter.atEnd(); ++inIter
623 : ) {
624 0 : mask(inIter.position()) = anyTrue(inIter.getMask());
625 : }
626 0 : _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 0 : }
635 :
636 0 : void ImageProfileFitterResults::_writeImages(
637 : const CoordinateSystem& xcsys,
638 : const Array<Bool>& mask, const String& yUnit
639 : ) const {
640 0 : map<String, String> mymap;
641 0 : map<String, String> unitmap;
642 0 : mymap["center"] = _centerName;
643 0 : mymap["centerErr"] = _centerErrName;
644 0 : mymap["fwhm"] = _fwhmName;
645 0 : mymap["fwhmErr"] = _fwhmErrName;
646 0 : mymap["amp"] = _ampName;
647 0 : mymap["ampErr"] = _ampErrName;
648 0 : mymap["integral"] = _integralName;
649 0 : mymap["integralErr"] = _integralErrName;
650 0 : unitmap["center"] = _xUnit;
651 0 : unitmap["centerErr"] = _xUnit;
652 0 : unitmap["fwhm"] = _xUnit;
653 0 : unitmap["fwhmErr"] = _xUnit;
654 0 : unitmap["amp"] = yUnit;
655 0 : unitmap["ampErr"] = yUnit;
656 0 : unitmap["integral"] = _xUnit + "." + yUnit;
657 0 : unitmap["integralErr"] = _xUnit + "." + yUnit;
658 0 : for (uInt i=0; i<_nGaussMultiplets+_nOthers; i++) {
659 0 : if (i == _gsPlane && _nGaussSinglets == 0) {
660 0 : continue;
661 : }
662 0 : else if (i == _lsPlane && _nLorentzSinglets == 0) {
663 0 : continue;
664 : }
665 0 : String id = _getTag(i);
666 0 : for (const auto& p : mymap) {
667 0 : String imagename = p.second;
668 : String suffix = i == _gsPlane
669 : ? ""
670 : : i == _lsPlane
671 : ? "_ls"
672 0 : : _nGaussMultiplets <= 1
673 : ? "_gm"
674 0 : : "_gm" + String::toString(i-_nOthers);
675 0 : imagename += suffix;
676 0 : if (! p.second.empty()) {
677 0 : _makeSolutionImages(
678 : imagename, xcsys,
679 0 : _results.asRecord(id).asArrayDouble(p.first),
680 0 : unitmap.find(p.first)->second, mask
681 : );
682 : }
683 : }
684 : }
685 0 : if (
686 0 : (_nPLPCoeffs > 0 && (! _plpName.empty() || ! _plpErrName.empty()))
687 0 : || (_nLTPCoeffs > 0 && (! _ltpName.empty() || ! _ltpErrName.empty()))
688 : ) {
689 0 : String type = _results.isDefined("plp") ? "plp" : "ltp";
690 0 : if (_nPLPCoeffs > 0 && ! _plpName.empty()) {
691 0 : _makeSolutionImages(
692 0 : _plpName, xcsys,
693 0 : _results.asRecord("plp").asArrayDouble("solution"),
694 : "", mask
695 : );
696 : }
697 0 : if (_nPLPCoeffs > 0 && ! _plpErrName.empty()) {
698 0 : _makeSolutionImages(
699 0 : _plpErrName, xcsys,
700 0 : _results.asRecord("plp").asArrayDouble("error"),
701 : "", mask
702 : );
703 : }
704 0 : if (_nLTPCoeffs > 0 && ! _ltpName.empty()) {
705 0 : _makeSolutionImages(
706 0 : _ltpName, xcsys,
707 0 : _results.asRecord("ltp").asArrayDouble("solution"),
708 : "", mask
709 : );
710 : }
711 0 : if (_nLTPCoeffs > 0 && ! _ltpErrName.empty()) {
712 0 : _makeSolutionImages(
713 0 : _ltpErrName, xcsys,
714 0 : _results.asRecord("ltp").asArrayDouble("error"),
715 : "", mask
716 : );
717 : }
718 : }
719 0 : }
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 0 : Double ImageProfileFitterResults::_fitAxisIncrement() const {
729 0 : if (_doVelocity) {
730 0 : Vector<Double> pixels(2);
731 0 : pixels[0] = 0;
732 0 : pixels[1] = 1;
733 0 : Vector<Double> velocities(2);
734 0 : _subImage->coordinates().spectralCoordinate().pixelToVelocity(
735 : velocities, pixels
736 : );
737 0 : return velocities[1] - velocities[0];
738 : }
739 : else {
740 0 : 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 0 : Double ImageProfileFitterResults::_centerWorld(
753 : const PCFSpectralElement& solution, const IPosition& imPos
754 : ) const {
755 0 : Vector<Double> pixel(imPos.size());
756 0 : for (uInt i=0; i<pixel.size(); i++) {
757 0 : pixel[i] = imPos[i];
758 : }
759 0 : Vector<Double> world(pixel.size());
760 : // in pixels here
761 0 : pixel[_fitAxis] = solution.getCenter();
762 0 : _subImage->coordinates().toWorld(world, pixel);
763 0 : if (_doVelocity) {
764 : Double velocity;
765 0 : _subImage->coordinates().spectralCoordinate().frequencyToVelocity(velocity, world(_fitAxis));
766 0 : return velocity;
767 : }
768 : else {
769 0 : return world(_fitAxis);
770 : }
771 : }
772 :
773 0 : void ImageProfileFitterResults::_resultsToLog() {
774 0 : if (! _logResults && _logfile.get() == 0) {
775 0 : return;
776 : }
777 0 : ostringstream summary;
778 0 : summary << "****** Fit performed at " << Time().toString() << "******" << endl << endl;
779 0 : summary << _summaryHeader;
780 0 : if (_nPLPCoeffs + _nLTPCoeffs == 0) {
781 0 : if (_goodAmpRange.size() == 2) {
782 0 : summary << " --- valid amplitude range: " << _goodAmpRange << endl;
783 : }
784 0 : if (_goodCenterRange.size() == 2) {
785 0 : summary << " --- valid center range: " << _goodCenterRange << endl;
786 : }
787 0 : if (_goodFWHMRange.size() == 2) {
788 0 : summary << " --- valid FWHM range: " << _goodFWHMRange << endl;
789 : }
790 0 : summary << " --- number of Gaussian singlets: " << _nGaussSinglets << endl;
791 0 : summary << " --- number of Gaussian multiplets: " << _nGaussMultiplets << endl;
792 0 : if (_nGaussMultiplets > 0) {
793 0 : for (uInt i=0; i<_nGaussMultiplets; i++) {
794 0 : Array<Double> amp = _results.asRecord("gm" + String::toString(i)).asArrayDouble(AMP);
795 0 : uInt n = amp.shape()[amp.ndim()-1];
796 0 : summary << " --- number of components in Gaussian multiplet "
797 0 : << i << ": " << n << endl;
798 : }
799 : }
800 0 : if (_polyOrder >= 0) {
801 0 : summary << " --- polynomial order: " << _polyOrder << endl;
802 : }
803 : else {
804 0 : summary << " --- no polynomial fit " << endl;
805 : }
806 : }
807 0 : if (_multiFit) {
808 0 : summary << " --- Multiple profiles fit, one per pixel over selected region" << endl;
809 : }
810 : else {
811 0 : summary << " --- One profile fit, averaged over several pixels if necessary" << endl;
812 : }
813 0 : if (_logResults) {
814 0 : *_log << LogIO::NORMAL << summary.str() << LogIO::POST;
815 : }
816 0 : _writeLogfile(summary.str(), False, False);
817 0 : IPosition inTileShape = _subImage->niceCursorShape();
818 0 : TiledLineStepper stepper (_subImage->shape(), inTileShape, _fitAxis);
819 0 : RO_MaskedLatticeIterator<Float> inIter(*_subImage, stepper);
820 0 : CoordinateSystem csysSub = _subImage->coordinates();
821 0 : Vector<Double> worldStart;
822 0 : ThrowIf(
823 : ! csysSub.toWorld(worldStart, inIter.position()),
824 : csysSub.errorMessage()
825 : );
826 0 : Vector<Double> pixStart;
827 0 : ThrowIf(
828 : ! _csysIm.toPixel(pixStart, worldStart),
829 : _csysIm.errorMessage()
830 : );
831 0 : if (_multiFit) {
832 0 : for (uInt i=0; i<pixStart.size(); i++) {
833 0 : pixStart[i] = (Int)std::floor( pixStart[i] + 0.5);
834 : }
835 : }
836 0 : Vector<Double> imPix(pixStart.size());
837 0 : Vector<Double> world;
838 0 : IPosition subimPos;
839 0 : SpectralList solutions;
840 0 : String axisUnit = _csysIm.worldAxisUnits()[_fitAxis];
841 0 : Int pixPrecision = _multiFit ? 0 : 3;
842 0 : _pixelPositions.resize( _fitters->size());
843 0 : uint index = 0;
844 0 : for (
845 0 : inIter.reset();
846 0 : ! inIter.atEnd();
847 0 : inIter++
848 : ) {
849 0 : subimPos = inIter.position();
850 0 : const std::shared_ptr<ProfileFitResults> fitter = (*_fitters)(subimPos);
851 0 : if (! fitter) {
852 0 : continue;
853 : }
854 0 : summary.str("");
855 0 : Int basePrecision = summary.precision(1);
856 0 : summary.precision(basePrecision);
857 0 : if (csysSub.toWorld(world, subimPos)) {
858 0 : summary << "Fit :" << endl;
859 0 : for (uInt i=0; i<world.size(); i++) {
860 0 : if ((Int)i != _fitAxis) {
861 0 : IPosition x(1, _axisTypes[i]);
862 0 : x.append(subimPos);
863 0 : if (_axisTypes[i] == LONGITUDE) {
864 0 : summary << " RA : "
865 0 : << _worldCoords(x) << endl;
866 : }
867 0 : else if (_axisTypes[i] == LATITUDE) {
868 0 : summary << " Dec : "
869 0 : << _worldCoords(x) << endl;
870 : }
871 0 : else if (_axisTypes[i] == FREQUENCY) {
872 0 : summary << " Freq : "
873 0 : << _worldCoords(x) << endl;
874 : }
875 0 : else if (_axisTypes[i] == POLARIZATION) {
876 0 : summary << " Stokes : " << _worldCoords(x) << endl;
877 : }
878 : }
879 : }
880 : }
881 : else {
882 0 : ThrowCc(csysSub.errorMessage());
883 : }
884 0 : for (uInt i=0; i<pixStart.size(); i++) {
885 0 : imPix[i] = pixStart[i] + subimPos[i];
886 : }
887 0 : _pixelPositions[index] = imPix;
888 0 : summary.setf(ios::fixed);
889 0 : summary << setprecision(pixPrecision) << " Pixel : [";
890 0 : for (uInt i=0; i<imPix.size(); i++) {
891 0 : if (i == (uInt)_fitAxis) {
892 0 : summary << " *";
893 : }
894 : else {
895 0 : summary << imPix[i];
896 : }
897 0 : if (i != imPix.size()-1) {
898 0 : summary << ", ";
899 : }
900 : }
901 0 : summary << "]" << setprecision(basePrecision) << endl;
902 0 : summary.unsetf(ios::fixed);
903 0 : Bool attempted = fitter->getList().nelements() > 0;
904 : summary << " Attempted : "
905 0 : << (attempted ? "YES" : "NO") << endl;
906 0 : if (attempted) {
907 0 : String converged = fitter->converged() ? "YES" : "NO";
908 0 : summary << " Converged : " << converged << endl;
909 0 : if (fitter->converged()) {
910 0 : summary << " Iterations : " << fitter->getNumberIterations() << endl;
911 0 : String valid = fitter->isValid() ? "YES" : "NO";
912 0 : summary << " Valid : " << valid << endl;
913 0 : if (fitter->isValid()) {
914 0 : solutions = fitter->getList();
915 0 : for (uInt i=0; i<solutions.nelements(); i++) {
916 0 : SpectralElement::Types type = solutions[i]->getType();
917 0 : if (solutions.nelements() > 1) {
918 0 : summary << " Results for component " << i << ":" << endl;
919 : }
920 0 : switch(type) {
921 0 : 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 0 : = dynamic_cast<const PCFSpectralElement*>(solutions[i]);
928 0 : summary << _pcfToString(
929 0 : pcf, _csysIm, world.copy(), subimPos
930 0 : );
931 : }
932 0 : break;
933 0 : case SpectralElement::GMULTIPLET:
934 : {
935 : const GaussianMultipletSpectralElement *gm
936 0 : = dynamic_cast<const GaussianMultipletSpectralElement*>(solutions[i]);
937 0 : summary << _gaussianMultipletToString(
938 0 : *gm, _csysIm, world.copy(), subimPos
939 0 : );
940 0 : 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 0 : case SpectralElement::POWERLOGPOLY:
950 : {
951 : const PowerLogPolynomialSpectralElement *p
952 0 : = dynamic_cast<const PowerLogPolynomialSpectralElement*>(solutions[i]);
953 0 : summary << _powerLogPolyToString(*p);
954 : }
955 0 : break;
956 0 : case SpectralElement::LOGTRANSPOLY:
957 : {
958 : const LogTransformedPolynomialSpectralElement *p
959 0 : = dynamic_cast<const LogTransformedPolynomialSpectralElement*>(solutions[i]);
960 0 : summary << _logTransPolyToString(*p);
961 : }
962 0 : break;
963 0 : default:
964 0 : ThrowCc("Logic Error: Unhandled spectral element type");
965 : break;
966 : }
967 : }
968 : }
969 : }
970 : }
971 0 : if (_logResults) {
972 0 : *_log << LogIO::NORMAL << summary.str() << endl << LogIO::POST;
973 : }
974 0 : _writeLogfile(summary.str(), False, False);
975 : }
976 0 : if (_logfile.get() != 0) {
977 0 : _logfile->close();
978 : }
979 : }
980 :
981 0 : String ImageProfileFitterResults::_elementToString(
982 : const Double value, const Double error,
983 : const String& unit, Bool isFixed
984 : ) const {
985 0 : Unit myUnit(unit);
986 0 : Vector<String> unitPrefix;
987 0 : String outUnit;
988 0 : Quantity qVal(value, unit);
989 0 : Quantity qErr(error, unit);
990 0 : 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 0 : else if (
1008 0 : unit.empty() || Quantity(1, myUnit).isConform(Quantity(1, "m/s"))
1009 0 : || isdigit(unit[0])
1010 : ) {
1011 : // do nothing
1012 : }
1013 : else {
1014 0 : Vector<String> unitPrefix(10);
1015 0 : unitPrefix[0] = "T";
1016 0 : unitPrefix[1] = "G";
1017 0 : unitPrefix[2] = "M";
1018 0 : unitPrefix[3] = "k";
1019 0 : unitPrefix[4] = "";
1020 0 : unitPrefix[5] = "m";
1021 0 : unitPrefix[6] = "u";
1022 0 : unitPrefix[7] = "n";
1023 0 : unitPrefix[8] = "p";
1024 0 : unitPrefix[9] = "f";
1025 :
1026 0 : for (auto prefix: unitPrefix) {
1027 0 : outUnit = prefix + unit;
1028 0 : if (fabs(qVal.getValue(outUnit)) > 1) {
1029 0 : qVal.convert(outUnit);
1030 0 : qErr.convert(outUnit);
1031 0 : break;
1032 : }
1033 : }
1034 : }
1035 0 : Vector<Double> valErr(2);
1036 0 : valErr[0] = qVal.getValue();
1037 0 : valErr[1] = qErr.getValue();
1038 :
1039 0 : uInt precision = precisionForValueErrorPairs(valErr, Vector<Double>());
1040 0 : ostringstream out;
1041 0 : out << std::fixed << setprecision(precision);
1042 0 : out << qVal.getValue();
1043 0 : if (isFixed && qErr.getValue() == 0) {
1044 0 : out << " (fixed)";
1045 : }
1046 : else {
1047 0 : out << " +/- " << qErr.getValue();
1048 : }
1049 0 : out << " " << qVal.getUnit();
1050 0 : return out.str();
1051 : }
1052 :
1053 0 : 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 0 : Vector<Double> myWorld = world;
1059 0 : String yUnit = _subImage->units().getName();
1060 0 : ostringstream summary;
1061 0 : Vector<Bool> fixed = pcf->fixed();
1062 0 : if (showTypeString) {
1063 0 : summary << indent << " Type : "
1064 0 : << SpectralElement::fromType(pcf->getType()) << endl;
1065 : }
1066 0 : summary << indent << " Peak : "
1067 0 : << _elementToString(
1068 0 : pcf->getAmpl(), pcf->getAmplErr(), yUnit, fixed[0]
1069 0 : ) << endl;
1070 0 : Double center = _centerWorld(
1071 : *pcf, imPos
1072 : );
1073 :
1074 0 : Double increment = fabs(_fitAxisIncrement());
1075 :
1076 0 : Double centerErr = pcf->getCenterErr() * increment;
1077 0 : Double fwhm = pcf->getFWHM() * increment;
1078 0 : Double fwhmErr = pcf->getFWHMErr() * increment;
1079 :
1080 0 : Double pCenter = 0;
1081 0 : Double pCenterErr = 0;
1082 0 : Double pFWHM = 0;
1083 0 : Double pFWHMErr = 0;
1084 0 : Int specCoordIndex = csys.findCoordinate(Coordinate::SPECTRAL);
1085 0 : Bool convertedCenterToPix = True;
1086 0 : Bool convertedFWHMToPix = True;
1087 0 : if (_doVelocity) {
1088 0 : if (csys.spectralCoordinate(specCoordIndex).velocityToPixel(pCenter, center)) {
1089 : Double nextVel;
1090 0 : csys.spectralCoordinate(specCoordIndex).pixelToVelocity(nextVel, pCenter+1);
1091 0 : Double velInc = fabs(center - nextVel);
1092 0 : pCenterErr = centerErr/velInc;
1093 0 : pFWHM = abs(fwhm/velInc);
1094 0 : 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 0 : summary << indent << " Center : "
1116 0 : << _elementToString(
1117 0 : center, centerErr, _xUnit, fixed[1]
1118 0 : ) << endl;
1119 0 : if (convertedCenterToPix) {
1120 0 : summary << indent << " "
1121 0 : << _elementToString(
1122 0 : pCenter, pCenterErr, "pixel", fixed[1]
1123 0 : ) << endl;
1124 : }
1125 : else {
1126 0 : summary << indent << " Could not convert world to pixel for center" << endl;
1127 : }
1128 0 : summary << indent << " FWHM : "
1129 0 : << _elementToString(
1130 0 : fwhm, fwhmErr, _xUnit, fixed[2]
1131 0 : )
1132 0 : << endl;
1133 0 : if (convertedFWHMToPix) {
1134 0 : summary << indent << " "
1135 0 : << _elementToString(pFWHM, pFWHMErr, "pixel", fixed[2])
1136 0 : << endl;
1137 : }
1138 : else {
1139 0 : summary << indent << " Could not convert FWHM to pixel" << endl;
1140 : }
1141 0 : Double integral = pcf->getIntegral()*increment;
1142 0 : Double integralErr = pcf->getIntegralErr()*increment;
1143 0 : String integUnit = (Quantity(1.0 ,yUnit)*Quantity(1.0, _xUnit)).getUnit();
1144 0 : summary << indent << " Integral : "
1145 0 : << _elementToString(integral, integralErr, integUnit, fixed[0] && fixed[2])
1146 0 : << endl;
1147 0 : if (fwhm/increment <= 3) {
1148 0 : summary << indent << "WARNING: The FWHM is only " << (fwhm/increment)
1149 : << " times the channel width. Be aware that instrumental channelization effects"
1150 0 : << " may be important." << endl;
1151 : }
1152 0 : return summary.str();
1153 : }
1154 :
1155 0 : String ImageProfileFitterResults::_gaussianMultipletToString(
1156 : const GaussianMultipletSpectralElement& gm,
1157 : const CoordinateSystem& csys,
1158 : const Vector<Double>& world, const IPosition& imPos
1159 : ) const {
1160 0 : Vector<GaussianSpectralElement> g(gm.getGaussians());
1161 0 : ostringstream summary;
1162 0 : summary << " Type : GAUSSIAN MULTIPLET" << endl;
1163 0 : for (uInt i=0; i<g.size(); i++) {
1164 0 : summary << " Results for subcomponent "
1165 0 : << i << ":" << endl;
1166 : summary
1167 0 : << _pcfToString(&g[i], csys, world, imPos, False, " ")
1168 0 : << endl;
1169 : }
1170 0 : 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 0 : String ImageProfileFitterResults::_powerLogPolyToString(
1243 : const PowerLogPolynomialSpectralElement& plp
1244 : ) const {
1245 0 : ostringstream summary;
1246 0 : summary << " Type : POWER LOGARITHMIC POLYNOMIAL" << endl;
1247 0 : summary << " Function : c0*(x/D)**(c1";
1248 0 : uInt nElements = plp.get().size();
1249 0 : 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 0 : summary << ")" << endl;
1256 0 : summary << " D : " << _plpDivisor << endl;
1257 0 : Vector<Double> parms = plp.get();
1258 0 : Vector<Double> errs = plp.getError();
1259 0 : Vector<Bool> fixed = plp.fixed();
1260 :
1261 0 : for (uInt j=0; j<parms.size(); j++) {
1262 0 : summary << " c" << j << " : "
1263 0 : << _elementToString(parms[j], errs[j], "", fixed[j]) << endl;
1264 : }
1265 0 : return summary.str();
1266 : }
1267 :
1268 0 : String ImageProfileFitterResults::_logTransPolyToString(
1269 : const LogTransformedPolynomialSpectralElement& ltp
1270 : ) const {
1271 0 : ostringstream summary;
1272 0 : summary << " Type : LOGARITHMIC TRANSFORMED POLYNOMIAL" << endl;
1273 0 : summary << " Function : ln(y) = c0";
1274 0 : uInt nElements = ltp.get().size();
1275 0 : for (uInt i=1; i<nElements; i++) {
1276 0 : summary << " + c" << i << "*ln(x/D)";
1277 0 : if (i > 2) {
1278 0 : summary << "**" << i;
1279 : }
1280 : }
1281 0 : summary << ")" << endl;
1282 0 : summary << " D : " << _plpDivisor << endl;
1283 0 : Vector<Double> parms = ltp.get();
1284 0 : Vector<Double> errs = ltp.getError();
1285 0 : Vector<Bool> fixed = ltp.fixed();
1286 :
1287 0 : for (uInt j=0; j<parms.size(); j++) {
1288 0 : summary << " c" << j << " : "
1289 0 : << _elementToString(parms[j], errs[j], "", fixed[j]) << endl;
1290 : }
1291 0 : return summary.str();
1292 : }
1293 :
1294 0 : 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 0 : auto valuesShape = values.shape();
1300 0 : Vector<Float> dataCopy(values.size());
1301 0 : Vector<Double>::const_iterator iter;
1302 : // isNaN(Array<Double>&) works, isNaN(Array<Float>&) gives spurious results
1303 0 : Array<Bool> nanInfMask = ! (isNaN(values) || isInf(values));
1304 0 : Vector<Float>::iterator jiter = dataCopy.begin();
1305 0 : for (iter=values.begin(); iter!=values.end(); ++iter, ++jiter) {
1306 0 : *jiter = (Float)*iter;
1307 : }
1308 0 : auto dc = dataCopy.reform(valuesShape);
1309 0 : uInt nImages = valuesShape.last();
1310 0 : auto imageShape = valuesShape.getFirst(valuesShape.size() - 1);
1311 0 : IPosition start(values.ndim(), 0);
1312 0 : IPosition end = valuesShape - 1;
1313 0 : end[end.size() - 1] = 0;
1314 0 : for (uInt i=0; i<nImages; ++i) {
1315 : auto myname = nImages == 0
1316 : ? name
1317 0 : : name + "_" + String::toString(i);
1318 0 : PagedImage<Float> image(imageShape, csys, myname);
1319 0 : start[start.size() - 1] = i;
1320 0 : end[end.size() - 1] = i;
1321 0 : auto imageVals = dc(start, end);
1322 0 : 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 0 : Array<Bool> nanInfMask = ! (isNaN(doubleValues) || isInf(imageVals));
1326 : // remove the last axis
1327 0 : imageVals.removeDegenerate(values.ndim() -1);
1328 0 : image.put(imageVals);
1329 0 : nanInfMask.removeDegenerate(values.ndim() -1);
1330 0 : Bool hasPixMask = ! allTrue(mask);
1331 0 : Bool hasNanMask = ! allTrue(nanInfMask);
1332 0 : if (hasNanMask || hasPixMask) {
1333 0 : Array<Bool> resMask(imageShape);
1334 : String maskName = image.makeUniqueRegionName(
1335 0 : String("mask"), 0
1336 0 : );
1337 0 : image.makeMask(maskName, True, True, False);
1338 0 : if (hasPixMask) {
1339 0 : resMask = mask;
1340 0 : if (hasNanMask) {
1341 0 : resMask = resMask && nanInfMask;
1342 : }
1343 : }
1344 : else {
1345 0 : resMask = nanInfMask;
1346 : }
1347 0 : (&image.pixelMask())->put(resMask);
1348 : }
1349 0 : image.setUnits(Unit(unit));
1350 : }
1351 0 : }
1352 :
1353 : }
|
