Line data Source code
1 :
2 : //# Copyright (C) 1998,1999,2000,2001,2003
3 : //# Associated Universities, Inc. Washington DC, USA.
4 : //#
5 : //# This program is free software; you can redistribute it and/or modify it
6 : //# under the terms of the GNU General Public License as published by the Free
7 : //# Software Foundation; either version 2 of the License, or (at your option)
8 : //# any later version.
9 : //#
10 : //# This program is distributed in the hope that it will be useful, but WITHOUT
11 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
13 : //# more details.
14 : //#
15 : //# You should have received a copy of the GNU General Public License along
16 : //# with this program; if not, write to the Free Software Foundation, Inc.,
17 : //# 675 Massachusetts Ave, Cambridge, MA 02139, USA.
18 : //#
19 : //# Correspondence concerning AIPS++ should be addressed as follows:
20 : //# Internet email: casa-feedback@nrao.edu.
21 : //# Postal address: AIPS++ Project Office
22 : //# National Radio Astronomy Observatory
23 : //# 520 Edgemont Road
24 : //# Charlottesville, VA 22903-2475 USA
25 : //#
26 : //# $Id: $
27 :
28 : #include <imageanalysis/ImageAnalysis/ImageProfileFitter.h>
29 :
30 : #include <casacore/casa/Quanta/MVAngle.h>
31 : #include <casacore/casa/Quanta/MVTime.h>
32 : #include <casacore/images/Images/ImageUtilities.h>
33 : #include <casacore/images/Images/PagedImage.h>
34 : #include <casacore/images/Images/TempImage.h>
35 : #include <casacore/scimath/Mathematics/Combinatorics.h>
36 : #include <casacore/casa/Arrays/ArrayLogical.h>
37 :
38 : #include <imageanalysis/ImageAnalysis/ProfileFitResults.h>
39 :
40 : #include <imageanalysis/ImageAnalysis/ImageCollapser.h>
41 : #include <imageanalysis/ImageAnalysis/SubImageFactory.h>
42 : #include <imageanalysis/IO/ProfileFitterEstimatesFileParser.h>
43 : #include <imageanalysis/IO/ImageProfileFitterResults.h>
44 :
45 : // debug
46 : #include <casacore/casa/OS/PrecTimer.h>
47 :
48 : using namespace casacore;
49 : namespace casa {
50 :
51 : const String ImageProfileFitter::_class = "ImageProfileFitter";
52 :
53 0 : ImageProfileFitter::ImageProfileFitter(
54 : const SPCIIF image, const String& region,
55 : const Record *const ®ionPtr, const String& box,
56 : const String& chans, const String& stokes,
57 : const String& mask, const Int axis,
58 : const uInt ngauss, Bool overwrite
59 0 : ) : ImageTask<Float>(
60 : image, region, regionPtr, box, chans, stokes,
61 : mask, "", False
62 : ),
63 0 : _residual(), _model(), _xUnit(), _centerName(),
64 0 : _centerErrName(), _fwhmName(), _fwhmErrName(),
65 0 : _ampName(), _ampErrName(), _integralName(),
66 0 : _integralErrName(), _plpName(), _plpErrName(), _sigmaName(),
67 0 : _abscissaDivisorForDisplay("1"), _multiFit(False),
68 0 : /*_deleteImageOnDestruct(False),*/ _logResults(True),
69 0 : _isSpectralIndex(False), _createResid(False), _overwrite(overwrite),
70 0 : _storeFits(True),
71 0 : _polyOrder(-1), _fitAxis(axis), _nGaussSinglets(ngauss),
72 0 : _nGaussMultiplets(0), _nLorentzSinglets(0), _nPLPCoeffs(0),
73 0 : _nLTPCoeffs(0), _minGoodPoints(1), _nProfiles(0), _nAttempted(0), _nSucceeded(0),
74 0 : _nConverged(0), _nValid(0), _results(Record()), _nonPolyEstimates(),
75 0 : _goodAmpRange(), _goodCenterRange(), _goodFWHMRange(),
76 0 : _sigma(), _abscissaDivisor(1.0), _residImage(), _goodPlanes() {
77 0 : *_getLog() << LogOrigin(_class, __func__);
78 0 : _construct();
79 0 : _finishConstruction();
80 0 : }
81 :
82 0 : ImageProfileFitter::ImageProfileFitter(
83 : const SPCIIF image, const String& region,
84 : const Record *const ®ionPtr, const String& box,
85 : const String& chans, const String& stokes,
86 : const String& mask, const Int axis,
87 : const String& estimatesFilename, Bool overwrite
88 0 : ) : ImageTask<Float>(
89 : image, region, regionPtr, box, chans, stokes,
90 : mask, "", False
91 : ),
92 0 : _residual(), _model(), _xUnit(), _centerName(),
93 0 : _centerErrName(), _fwhmName(), _fwhmErrName(),
94 0 : _ampName(), _ampErrName(), _integralName(),
95 0 : _integralErrName(), _plpName(), _plpErrName(), _sigmaName(),
96 0 : _abscissaDivisorForDisplay("1"), _multiFit(False),
97 0 : /*_deleteImageOnDestruct(False),*/ _logResults(True),
98 0 : _isSpectralIndex(False), _createResid(False), _overwrite(overwrite),
99 0 : _storeFits(True),
100 0 : _polyOrder(-1), _fitAxis(axis), _nGaussSinglets(0),
101 0 : _nGaussMultiplets(0), _nLorentzSinglets(0), _nPLPCoeffs(0),
102 0 : _nLTPCoeffs(0), _minGoodPoints(1), _nProfiles(0), _nAttempted(0), _nSucceeded(0),
103 0 : _nConverged(0), _nValid(0), _results(Record()), _nonPolyEstimates(),
104 0 : _goodAmpRange(), _goodCenterRange(), _goodFWHMRange(),
105 0 : _sigma(), _abscissaDivisor(1.0), _residImage(), _goodPlanes() {
106 0 : *_getLog() << LogOrigin(_class, __func__);
107 0 : ThrowIf(estimatesFilename.empty(), "Estimates filename cannot be empty");
108 0 : ProfileFitterEstimatesFileParser parser(estimatesFilename);
109 0 : _nonPolyEstimates = parser.getEstimates();
110 0 : _nGaussSinglets = _nonPolyEstimates.nelements();
111 0 : *_getLog() << LogIO::NORMAL << "Number of gaussian singlets to fit found to be "
112 : <<_nGaussSinglets << " in estimates file " << estimatesFilename
113 0 : << LogIO::POST;
114 0 : _construct();
115 0 : _finishConstruction();
116 0 : }
117 :
118 0 : ImageProfileFitter::ImageProfileFitter(
119 : const SPCIIF image, const String& region,
120 : const Record *const ®ionPtr, const String& box,
121 : const String& chans, const String& stokes,
122 : const String& mask, const Int axis,
123 : const SpectralList& spectralList, Bool overwrite
124 0 : ) : ImageTask<Float>(
125 : image, region, regionPtr, box, chans, stokes,
126 : mask, "", False
127 : ),
128 0 : _residual(), _model(), _xUnit(), _centerName(),
129 0 : _centerErrName(), _fwhmName(), _fwhmErrName(),
130 0 : _ampName(), _ampErrName(), _integralName(),
131 0 : _integralErrName(), _plpName(), _plpErrName(), _sigmaName(),
132 0 : _abscissaDivisorForDisplay("1"), _multiFit(False),
133 0 : /* _deleteImageOnDestruct(False),*/ _logResults(True),
134 0 : _isSpectralIndex(False), _createResid(False), _overwrite(overwrite),
135 0 : _storeFits(True),
136 0 : _polyOrder(-1), _fitAxis(axis), _nGaussSinglets(0),
137 0 : _nGaussMultiplets(0), _nLorentzSinglets(0), _nPLPCoeffs(0),
138 0 : _nLTPCoeffs(0), _minGoodPoints(1), _nProfiles(0), _nAttempted(0), _nSucceeded(0),
139 0 : _nConverged(0), _nValid(0), _results(Record()), _nonPolyEstimates(),
140 0 : _goodAmpRange(), _goodCenterRange(), _goodFWHMRange(),
141 0 : _sigma(), _abscissaDivisor(1.0), _residImage(), _goodPlanes() {
142 0 : *_getLog() << LogOrigin(_class, __func__);
143 0 : ThrowIf(
144 : spectralList.nelements() == 0,
145 : "spectralList cannot be empty"
146 : );
147 0 : _nonPolyEstimates = spectralList;
148 0 : _nGaussSinglets = 0;
149 0 : _nGaussMultiplets = 0;
150 0 : for (uInt i=0; i<_nonPolyEstimates.nelements(); i++) {
151 0 : SpectralElement::Types myType = _nonPolyEstimates[i]->getType();
152 0 : switch(myType) {
153 0 : case SpectralElement::GAUSSIAN:
154 0 : _nGaussSinglets++;
155 0 : break;
156 0 : case SpectralElement::GMULTIPLET:
157 0 : _nGaussMultiplets++;
158 0 : break;
159 0 : case SpectralElement::LORENTZIAN:
160 0 : _nLorentzSinglets++;
161 0 : break;
162 0 : case SpectralElement::POWERLOGPOLY:
163 0 : ThrowIf(
164 : _nonPolyEstimates.nelements() > 1 || _polyOrder > 0,
165 : "Only a single power logarithmic polynomial may be fit "
166 : "and it cannot be fit simultaneously with other functions"
167 : );
168 0 : _nPLPCoeffs = _nonPolyEstimates[i]->get().size();
169 0 : break;
170 0 : case SpectralElement::LOGTRANSPOLY:
171 0 : ThrowIf(
172 : _nonPolyEstimates.nelements() > 1 || _polyOrder > 0,
173 : "Only a single transformed logarithmic polynomial may "
174 : "be fit and it cannot be fit simultaneously with other functions"
175 : );
176 0 : _nLTPCoeffs = _nonPolyEstimates[i]->get().size();
177 0 : break;
178 0 : default:
179 0 : ThrowCc(
180 : "Logic error: Only Gaussian singlets, "
181 : "Gaussian multiplets, and Lorentzian singlets, or a single power "
182 : "logarithmic polynomial, or a single log transformed polynomial are "
183 : "permitted in the spectralList input parameter"
184 : );
185 : break;
186 : }
187 0 : if (_nPLPCoeffs > 0) {
188 0 : *_getLog() << LogIO::NORMAL << "Will fit a single power logarithmic polynomial "
189 0 : << " from provided spectral element list" << LogIO::POST;
190 : }
191 0 : else if (_nLTPCoeffs > 0) {
192 0 : *_getLog() << LogIO::NORMAL << "Will fit a single logarithmic transformed polynomial "
193 0 : << " from provided spectral element list" << LogIO::POST;
194 : }
195 : else {
196 0 : if (_nGaussSinglets > 0) {
197 0 : *_getLog() << LogIO::NORMAL << "Number of Gaussian singlets to fit found to be "
198 : << _nGaussSinglets << " from provided spectral element list"
199 0 : << LogIO::POST;
200 : }
201 0 : if (_nGaussMultiplets > 0) {
202 0 : *_getLog() << LogIO::NORMAL << "Number of Gaussian multiplets to fit found to be "
203 : << _nGaussMultiplets << " from provided spectral element list"
204 0 : << LogIO::POST;
205 : }
206 0 : if (_nLorentzSinglets > 0) {
207 0 : *_getLog() << LogIO::NORMAL << "Number of lorentzian singlets to fit found to be "
208 : << _nLorentzSinglets << " from provided spectral element list"
209 0 : << LogIO::POST;
210 : }
211 : }
212 : }
213 0 : _construct();
214 0 : _finishConstruction();
215 0 : }
216 :
217 0 : ImageProfileFitter::~ImageProfileFitter() {}
218 :
219 0 : Record ImageProfileFitter::fit(Bool doDetailedResults) {
220 : // do this check here rather than at construction because _polyOrder can be set
221 : // after construction but before fit() is called
222 0 : _checkNGaussAndPolyOrder();
223 0 : LogOrigin logOrigin(_class, __func__);
224 0 : *_getLog() << logOrigin;
225 0 : std::unique_ptr<ImageInterface<Float> > originalSigma;
226 : {
227 0 : _subImage = SubImageFactory<Float>::createSubImageRO(
228 0 : *_getImage(), *_getRegion(), _getMask(), 0,
229 0 : AxesSpecifier(), _getStretch()
230 0 : );
231 0 : uInt nUnknowns = _nUnknowns();
232 0 : ThrowIf(
233 : nUnknowns >= _subImage->shape()[_fitAxis],
234 : "There are not enough points ("
235 : + String::toString(_subImage->shape()[_fitAxis])
236 : + ") along the fit axis to fit " + String::toString(nUnknowns)
237 : + " unknowns"
238 : );
239 0 : if (_sigma.get()) {
240 0 : if (! _sigmaName.empty()) {
241 0 : originalSigma.reset(_sigma->cloneII());
242 : }
243 : std::shared_ptr<const SubImage<Float> > sigmaSubImage = SubImageFactory<Float>::createSubImageRO(
244 0 : *_sigma, *_getRegion(), _getMask(), 0, AxesSpecifier(), _getStretch()
245 0 : );
246 0 : _sigma.reset(
247 : new TempImage<Float>(
248 0 : sigmaSubImage->shape(), sigmaSubImage->coordinates()
249 0 : )
250 : );
251 0 : _sigma->put(sigmaSubImage->get());
252 0 : }
253 : }
254 0 : *_getLog() << logOrigin;
255 0 : _storeFits = doDetailedResults || ! _centerName.empty()
256 0 : || ! _centerErrName.empty() || ! _fwhmName.empty()
257 0 : || ! _fwhmErrName.empty() || ! _ampName.empty()
258 0 : || ! _ampErrName.empty() || ! _integralName.empty()
259 0 : || ! _integralErrName.empty()
260 0 : || ! _plpName.empty() || ! _plpErrName.empty()
261 0 : || ! _ltpName.empty() || ! _ltpErrName.empty();
262 : try {
263 0 : if (! _multiFit) {
264 : ImageCollapser<Float> collapser(
265 0 : _subImage, IPosition(1, _fitAxis), True,
266 : ImageCollapserData::MEAN, "", True
267 0 : );
268 0 : SPIIF x = collapser.collapse();
269 : // _subImage needs to be a SubImage<Float> object
270 0 : _subImage = SubImageFactory<Float>::createSubImageRO(
271 0 : *x, Record(), "", _getLog().get(),
272 0 : AxesSpecifier(), False
273 0 : );
274 0 : if (_sigma.get()) {
275 0 : Array<Bool> sigmaMask = _sigma->get() != Array<Float>(_sigma->shape(), 0.0f);
276 0 : if (anyTrue(! sigmaMask)) {
277 0 : if (_sigma->hasPixelMask()) {
278 0 : sigmaMask = sigmaMask && _sigma->pixelMask().get();
279 : }
280 : else {
281 0 : _sigma->makeMask("sigmamask", True, True, False);
282 : }
283 0 : _sigma->pixelMask().put(sigmaMask);
284 : }
285 : ImageCollapser<Float> collapsedSigma(
286 0 : _sigma, IPosition(1, _fitAxis), True,
287 : ImageCollapserData::MEAN, "", True
288 0 : );
289 0 : SPIIF collapsed = collapsedSigma.collapse();
290 0 : std::shared_ptr<TempImage<Float> >ctmp = std::dynamic_pointer_cast<TempImage<Float> >(collapsed);
291 0 : ThrowIf(! ctmp, "Dynamic cast failed");
292 0 : _sigma = ctmp;
293 0 : }
294 0 : }
295 0 : _fitallprofiles();
296 0 : *_getLog() << logOrigin;
297 : }
298 0 : catch (const AipsError& x) {
299 0 : ThrowCc("Exception during fit: " + x.getMesg());
300 0 : }
301 : ImageProfileFitterResults resultHandler(
302 0 : _getLog(), _getImage()->coordinates(), &_fitters,
303 0 : _nonPolyEstimates, _subImage, _polyOrder,
304 : _fitAxis, _nGaussSinglets, _nGaussMultiplets,
305 0 : _nLorentzSinglets, _nPLPCoeffs, _nLTPCoeffs, _logResults,
306 0 : _multiFit, _getLogFile(), _xUnit, _summaryHeader()
307 0 : );
308 0 : addHistory(
309 : logOrigin,
310 0 : resultHandler.logSummary(
311 : _nProfiles, _nAttempted, _nSucceeded, _nConverged, _nValid
312 : )
313 : );
314 0 : if (_nPLPCoeffs > 0) {
315 0 : resultHandler.setPLPName(_plpName);
316 0 : resultHandler.setPLPErrName(_plpErrName);
317 0 : resultHandler.setPLPDivisor(_abscissaDivisorForDisplay);
318 : }
319 0 : else if (_nLTPCoeffs > 0) {
320 0 : resultHandler.setLTPName(_ltpName);
321 0 : resultHandler.setLTPErrName(_ltpErrName);
322 0 : resultHandler.setPLPDivisor(_abscissaDivisorForDisplay);
323 : }
324 0 : else if (_nGaussSinglets + _nGaussMultiplets + _nLorentzSinglets > 0) {
325 0 : resultHandler.setAmpName(_ampName);
326 0 : resultHandler.setAmpErrName(_ampErrName);
327 0 : resultHandler.setCenterName(_centerName);
328 0 : resultHandler.setCenterErrName(_centerErrName);
329 0 : resultHandler.setFWHMName(_fwhmName);
330 0 : resultHandler.setFWHMErrName(_fwhmErrName);
331 0 : resultHandler.setIntegralName(_integralName);
332 0 : resultHandler.setIntegralErrName(_integralErrName);
333 : }
334 0 : if (doDetailedResults) {
335 0 : resultHandler.createResults();
336 0 : _results = resultHandler.getResults();
337 : }
338 0 : resultHandler.writeImages(_nConverged > 0);
339 0 : if (_modelImage) {
340 0 : _modelImage = _prepareOutputImage(*_modelImage, _model, _overwrite, True);
341 : }
342 0 : if (_residImage) {
343 0 : _residImage = _prepareOutputImage(*_residImage, _residual, _overwrite, True);
344 : }
345 0 : if (originalSigma && ! _sigmaName.empty()) {
346 0 : _prepareOutputImage(*originalSigma, _sigmaName, True, True);
347 : }
348 0 : return _results;
349 0 : }
350 :
351 0 : uInt ImageProfileFitter::_nUnknowns() const {
352 0 : uInt n = 0;
353 0 : if (_polyOrder >= 0) {
354 0 : n += _polyOrder + 1;
355 : }
356 0 : if (_nGaussSinglets > 0) {
357 0 : n += 3*_nGaussSinglets;
358 : }
359 0 : uInt nel = _nonPolyEstimates.nelements();
360 0 : if (n == 0) {
361 0 : return n;
362 : }
363 0 : for (uInt i=0; i<nel; ++i) {
364 0 : const SpectralElement *const x = _nonPolyEstimates[i];
365 0 : Vector<Bool> fixed = x->fixed();
366 0 : Vector<Bool>::const_iterator iter = fixed.begin();
367 0 : Vector<Bool>::const_iterator end = fixed.end();
368 0 : while (iter != end) {
369 0 : if (*iter) {
370 0 : --n;
371 0 : if (n == 0) {
372 0 : return n;
373 : }
374 : }
375 0 : ++iter;
376 : }
377 0 : }
378 0 : return n;
379 : }
380 :
381 :
382 0 : void ImageProfileFitter::setPolyOrder(Int p) {
383 0 : ThrowIf(p < 0,"A polynomial cannot have a negative order");
384 0 : ThrowIf(
385 : _nPLPCoeffs > 0,
386 : "Cannot simultaneously fit a polynomial and "
387 : "a power logarithmic polynomial."
388 : );
389 0 : ThrowIf(
390 : _nLTPCoeffs > 0,
391 : "Cannot simultaneously fit a polynomial and "
392 : "a logarithmic transformed polynomial"
393 : );
394 0 : _polyOrder = p;
395 0 : }
396 :
397 0 : void ImageProfileFitter::setGoodAmpRange(const Double minv, const Double maxv) {
398 0 : _goodAmpRange.reset(
399 0 : new std::pair<Double, Double>(min(minv, maxv), max(minv, maxv))
400 : );
401 0 : }
402 :
403 0 : void ImageProfileFitter::setGoodCenterRange(const Double minv, const Double maxv) {
404 0 : _goodCenterRange.reset(
405 0 : new std::pair<Double, Double>(min(minv, maxv), max(minv, maxv))
406 : );
407 0 : }
408 :
409 0 : void ImageProfileFitter::setGoodFWHMRange(const Double minv, const Double maxv) {
410 0 : _goodFWHMRange.reset(
411 0 : new std::pair<Double, Double>(min(minv, maxv), max(minv, maxv))
412 : );
413 0 : }
414 :
415 0 : void ImageProfileFitter::setSigma(const Array<Float>& sigma) {
416 0 : std::unique_ptr<TempImage<Float> > temp;
417 0 : if (sigma.ndim() == _getImage()->ndim()) {
418 0 : temp.reset(new TempImage<Float>(
419 0 : sigma.shape(), _getImage()->coordinates())
420 : );
421 : }
422 0 : else if (sigma.ndim() == 1) {
423 0 : SpectralCoordinate sp;
424 0 : CoordinateSystem csys;
425 0 : csys.addCoordinate(sp);
426 0 : temp.reset(new TempImage<Float>(sigma.shape(), csys));
427 0 : }
428 0 : temp->put(sigma);
429 0 : setSigma(temp.get());
430 0 : }
431 :
432 0 : void ImageProfileFitter::setSigma(const ImageInterface<Float>* const &sigma) {
433 0 : if (anyTrue(sigma->get() < Array<Float>(sigma->shape(), 0.0f))) {
434 0 : *_getLog() << "All sigma values must be non-negative" << LogIO::EXCEPTION;
435 : }
436 0 : Float mymax = fabs(max(sigma->get()));
437 0 : if (
438 0 : sigma->ndim() == _getImage()->ndim()
439 0 : && sigma->shape() == _getImage()->shape()
440 : ) {
441 0 : SPIIF clone(sigma->cloneII());
442 0 : _sigma = std::dynamic_pointer_cast<TempImage<Float> >(clone);
443 0 : if (! _sigma) {
444 : SPIIF x = SubImageFactory<Float>::createImage(
445 0 : *sigma, "", Record(), "", False, False ,True, False
446 0 : );
447 0 : if (x) {
448 0 : _sigma = std::dynamic_pointer_cast<TempImage<Float> >(x);
449 0 : ThrowIf(
450 : ! _sigma,
451 : "Unable to create temporary weights image"
452 : );
453 : }
454 0 : }
455 0 : if (mymax != 1) {
456 0 : _sigma->put(_sigma->get()/mymax);
457 : }
458 0 : }
459 0 : else if (
460 0 : sigma->ndim() == _getImage()->ndim()
461 0 : || sigma->ndim() == 1
462 : ) {
463 0 : if (sigma->ndim() == _getImage()->ndim()) {
464 0 : IPosition expShape(_getImage()->ndim(), 1);
465 0 : expShape[_fitAxis] = _getImage()->shape()[_fitAxis];
466 0 : ThrowIf(
467 : sigma->shape() != expShape,
468 : "If the shape of the standard deviation image differs "
469 : "from the shape of the input image, the shape of the "
470 : "standard deviation image must be " + expShape.toString()
471 : );
472 0 : }
473 0 : else if (sigma->ndim() == 1) {
474 0 : ThrowIf(
475 : sigma->shape()[0] != _getImage()->shape()[_fitAxis],
476 : "A one dimensional standard deviation spectrum must have the same "
477 : "number of pixels as the input image has along its axis to be fit"
478 : );
479 : }
480 0 : IPosition dataToInsertShape(_getImage()->ndim(), 1);
481 0 : dataToInsertShape[_fitAxis] = _getImage()->shape()[_fitAxis];
482 0 : _sigma.reset(new TempImage<Float>(_getImage()->shape(), _getImage()->coordinates()));
483 0 : Array<Float> dataToInsert(IPosition(1, _getImage()->shape()[_fitAxis]));
484 0 : dataToInsert = sigma->get(sigma->ndim() == _getImage()->ndim());
485 : // normalize
486 0 : if (mymax != 1) {
487 0 : dataToInsert /= mymax;
488 : }
489 0 : Array<Float> sigmaData = _sigma->get();
490 0 : ArrayIterator<Float> iter(sigmaData, IPosition(1, _fitAxis), True);
491 0 : while(! iter.pastEnd()) {
492 0 : iter.array() = dataToInsert;
493 0 : iter.next();
494 : }
495 0 : _sigma->put(sigmaData);
496 0 : }
497 : else {
498 0 : ThrowCc("Illegal shape of standard deviation image");
499 : }
500 0 : if (! _sigma->coordinates().near(_getImage()->coordinates())) {
501 0 : _sigma->setCoordinateInfo(_getImage()->coordinates());
502 : }
503 0 : }
504 :
505 0 : Record ImageProfileFitter::getResults() const {
506 0 : return _results;
507 : }
508 :
509 0 : void ImageProfileFitter::setAbscissaDivisor(Double d) {
510 0 : if (! _isSpectralIndex) {
511 0 : *_getLog() << LogOrigin(_class, __func__);
512 0 : *_getLog() << LogIO::WARN << "This object is not configured to fit a "
513 : << "spectral index function, and so setting the abscissa divisor "
514 0 : << "will have no effect in the fitting process." << LogIO::POST;
515 : }
516 0 : _abscissaDivisor = d;
517 0 : _abscissaDivisorForDisplay = String::toString(d)
518 0 : + _getImage()->coordinates().worldAxisUnits()[_fitAxis];
519 0 : }
520 :
521 0 : void ImageProfileFitter::setAbscissaDivisor(const Quantity& q) {
522 0 : String fitAxisUnit = _getImage()->coordinates().worldAxisUnits()[_fitAxis];
523 0 : ThrowIf(
524 : ! q.isConform(fitAxisUnit),
525 : "Abscissa divisor unit " + q.getUnit()
526 : + " is not consistent with fit axis unit"
527 : );
528 0 : if (! _isSpectralIndex) {
529 0 : *_getLog() << LogOrigin(_class, __func__);
530 0 : *_getLog() << LogIO::WARN << "This object is not configured to fit a spectral index function "
531 : << "and so setting the abscissa divisor will have no effect in the fitting process."
532 0 : << LogIO::POST;
533 : }
534 0 : _abscissaDivisor = q.getValue(fitAxisUnit);
535 0 : _abscissaDivisorForDisplay = String::toString(q);
536 0 : }
537 :
538 0 : std::vector<OutputDestinationChecker::OutputStruct> ImageProfileFitter::_getOutputStruct() {
539 0 : std::vector<OutputDestinationChecker::OutputStruct> outputs;
540 0 : if (! _model.empty()) {
541 0 : OutputDestinationChecker::OutputStruct modelImage;
542 0 : modelImage.label = "model image";
543 0 : modelImage.outputFile = &_model;
544 0 : modelImage.required = True;
545 0 : modelImage.replaceable = _overwrite;
546 0 : outputs.push_back(modelImage);
547 0 : }
548 0 : if (! _residual.empty()) {
549 0 : OutputDestinationChecker::OutputStruct residImage;
550 0 : residImage.label = "residual image";
551 0 : residImage.outputFile = &_residual;
552 0 : residImage.required = True;
553 0 : residImage.replaceable = _overwrite;
554 0 : outputs.push_back(residImage);
555 0 : }
556 0 : return outputs;
557 0 : }
558 :
559 0 : void ImageProfileFitter::_checkNGaussAndPolyOrder() const {
560 0 : ThrowIf(
561 : _polyOrder < 0
562 : && (
563 : _nGaussSinglets + _nGaussMultiplets
564 : + _nLorentzSinglets
565 : ) == 0 && ! _isSpectralIndex,
566 : "Number of non-polynomials is 0 and polynomial order is less than zero. "
567 : "According to these inputs there is nothing to fit."
568 : );
569 0 : }
570 :
571 0 : void ImageProfileFitter::_finishConstruction() {
572 0 : LogOrigin logOrigin(_class, __func__);
573 0 : _isSpectralIndex = _nPLPCoeffs + _nLTPCoeffs > 0;
574 0 : ThrowIf(
575 : _fitAxis >= (Int)_getImage()->ndim(),
576 : "Specified fit axis " + String::toString(_fitAxis)
577 : + " must be less than the number of image axes ("
578 : + String::toString(_getImage()->ndim()) + ")"
579 : )
580 0 : if (_fitAxis < 0) {
581 0 : if (! _getImage()->coordinates().hasSpectralAxis()) {
582 0 : _fitAxis = 0;
583 0 : *_getLog() << LogIO::WARN << "No spectral coordinate found in image, "
584 0 : << "using axis 0 as fit axis" << LogIO::POST;
585 : }
586 : else {
587 0 : _fitAxis = _getImage()->coordinates().spectralAxisNumber();
588 0 : *_getLog() << LogIO::NORMAL << "Using spectral axis (axis " << _fitAxis
589 0 : << ") as fit axis" << LogIO::POST;
590 : }
591 : }
592 : //this->_setSupportsLogfile(true);
593 0 : }
594 : /*
595 : Bool ImageProfileFitter::_inVelocitySpace() const {
596 : return _fitAxis == _subImage->coordinates().spectralAxisNumber()
597 : && Quantity(1, _xUnit).isConform("m/s")
598 : && _subImage->coordinates().spectralCoordinate().restFrequency() > 0;
599 : }
600 : */
601 :
602 0 : Double ImageProfileFitter::getWorldValue(
603 : double pixelVal, const IPosition& imPos, const String& units,
604 : bool velocity, bool wavelength, int tabularIndex, MFrequency::Types freqType ) const {
605 0 : Vector<Double> pixel(imPos.size());
606 0 : for (uInt i=0; i<pixel.size(); i++) {
607 0 : pixel[i] = imPos[i];
608 : }
609 0 : Vector<Double> world(pixel.size());
610 : // in pixels here
611 0 : pixel[_fitAxis] = pixelVal;
612 0 : _subImage->coordinates().toWorld(world, pixel);
613 0 : SpectralCoordinate spectCoord;
614 : //Use a tabular index for conversion if one has been specified.
615 0 : if ( tabularIndex >= 0 ) {
616 0 : const CoordinateSystem& cSys = _subImage->coordinates();
617 0 : TabularCoordinate tabCoord = cSys.tabularCoordinate( tabularIndex );
618 0 : Vector<Double> worlds = tabCoord.worldValues();
619 0 : spectCoord = SpectralCoordinate( freqType, worlds );
620 0 : }
621 : //Use the default spectral axis for conversion
622 : else {
623 0 : spectCoord = _subImage->coordinates().spectralCoordinate();
624 : }
625 : Double convertedVal;
626 0 : if (velocity) {
627 0 : spectCoord.setVelocity( units );
628 0 : spectCoord.frequencyToVelocity(convertedVal, world(_fitAxis));
629 : }
630 0 : else if ( wavelength ) {
631 0 : spectCoord.setWavelengthUnit( units );
632 0 : Vector<Double> worldVal(1);
633 0 : worldVal[0] = world(_fitAxis);
634 0 : Vector<Double> waveVal(1);
635 0 : spectCoord.frequencyToWavelength(waveVal, worldVal);
636 0 : convertedVal = waveVal[0];
637 0 : }
638 : else {
639 0 : convertedVal = world(_fitAxis);
640 : }
641 0 : return convertedVal;
642 0 : }
643 :
644 0 : void ImageProfileFitter::_fitallprofiles() {
645 0 : *_getLog() << LogOrigin(_class, __func__);
646 : // Create output images with a mask. Make them TempImages to start
647 : // in attempt to improve IO performance, and write them out if necessary
648 : // at the end
649 0 : if (
650 0 : ! _model.empty()
651 0 : && ! (
652 0 : _modelImage = SubImageFactory<Float>::createImage(
653 0 : *_subImage, "", Record(), "",
654 0 : False, _overwrite ,False, False, True
655 0 : )
656 0 : )
657 : ) {
658 0 : *_getLog() << LogIO::WARN << "Failed to create model image" << LogIO::POST;
659 : }
660 0 : if (
661 0 : (! _residual.empty() || _createResid)
662 0 : && ! (
663 0 : _residImage = SubImageFactory<Float>::createImage(
664 0 : *_subImage, "", Record(), "",
665 0 : False, _overwrite ,False, False, True
666 0 : )
667 0 : )
668 : ) {
669 0 : *_getLog() << LogIO::WARN << "Failed to create residual image" << LogIO::POST;
670 : }
671 0 : Bool showProgress = True;
672 0 : _fitProfiles(showProgress);
673 0 : }
674 :
675 0 : void ImageProfileFitter::_fitProfiles(Bool showProgress) {
676 0 : IPosition inShape = _subImage->shape();
677 0 : if (_modelImage) {
678 0 : AlwaysAssert(inShape.isEqual(_modelImage->shape()), AipsError);
679 : }
680 0 : if (_residImage) {
681 0 : AlwaysAssert(inShape.isEqual(_residImage->shape()), AipsError);
682 : }
683 0 : const auto nDim = _subImage->ndim();
684 0 : IPosition sliceShape(nDim, 1);
685 0 : sliceShape(_fitAxis) = inShape(_fitAxis);
686 0 : Array<Float> resultData(sliceShape);
687 0 : Array<Bool> resultMask(sliceShape);
688 0 : String doppler = "";
689 0 : auto csys = _subImage->coordinates();
690 0 : _xUnit = csys.spectralCoordinate().worldAxisUnits()[0];
691 0 : if (
692 0 : ! _isSpectralIndex && _fitAxis == _subImage->coordinates().spectralAxisNumber()
693 0 : && _subImage->coordinates().spectralCoordinate().restFrequency() > 0
694 : ) {
695 0 : SpectralCoordinate specCoord = csys.spectralCoordinate();
696 0 : _xUnit = specCoord.velocityUnit();
697 : doppler = MDoppler::showType(
698 : specCoord.velocityDoppler()
699 0 : );
700 0 : }
701 0 : String errMsg;
702 : ImageFit1D<Float>::AbcissaType abcissaType;
703 0 : auto abscissaUnits = _isSpectralIndex ? "native" : "pix";
704 0 : ThrowIf(
705 : ! ImageFit1D<Float>::setAbcissaState(
706 : errMsg, abcissaType, csys, abscissaUnits, doppler, _fitAxis
707 : ), errMsg
708 : );
709 0 : IPosition fitterShape = inShape;
710 0 : fitterShape[_fitAxis] = 1;
711 0 : if (_storeFits) {
712 0 : _fitters.resize(fitterShape);
713 : }
714 0 : _nProfiles = fitterShape.product();
715 0 : std::shared_ptr<ProgressMeter> pProgressMeter;
716 0 : if (showProgress) {
717 0 : ostringstream oss;
718 0 : oss << "Fit profiles on axis " << _fitAxis+1;
719 0 : pProgressMeter.reset(
720 0 : new ProgressMeter(0, _nProfiles, oss.str())
721 : );
722 0 : }
723 0 : SPCIIF fitData = _subImage;
724 0 : std::set<uInt> myGoodPlanes;
725 0 : if (! _goodPlanes.empty()) {
726 0 : IPosition origin(_subImage->ndim(), 0);
727 0 : Vector<Double> world(_subImage->ndim(), 0);
728 0 : csys.toWorld(world, origin);
729 0 : const CoordinateSystem imcsys = _getImage()->coordinates();
730 0 : Int imageOff = Int(imcsys.toPixel(world)[_fitAxis] + 0.5);
731 0 : AlwaysAssert(imageOff >= 0, AipsError);
732 0 : std::vector<Int> goodPlanes(_goodPlanes.begin(), _goodPlanes.end());
733 0 : if (imageOff > 0) {
734 0 : goodPlanes = std::vector<Int>(_goodPlanes.size());
735 0 : std::transform(
736 : _goodPlanes.begin(), _goodPlanes.end(), goodPlanes.begin(),
737 0 : bind2nd(minus<Int>(), imageOff)
738 : );
739 : }
740 0 : std::vector<Int>::iterator iter = goodPlanes.begin();
741 0 : while (iter != goodPlanes.end() && *iter < 0) {
742 0 : goodPlanes.erase(iter);
743 0 : iter = goodPlanes.begin();
744 : }
745 0 : myGoodPlanes = std::set<uInt>(goodPlanes.begin(), goodPlanes.end());
746 0 : }
747 0 : Bool checkMinPts = fitData->isMasked();
748 0 : Array<Bool> fitMask;
749 0 : if (checkMinPts) {
750 : fitMask = (
751 0 : partialNTrue(fitData->getMask(False), IPosition(1, _fitAxis))
752 0 : >= (long unsigned int) _minGoodPoints
753 0 : );
754 0 : IPosition oldShape = fitMask.shape();
755 0 : IPosition newShape(fitMask.ndim() + 1);
756 0 : uInt oldIndex = 0;
757 0 : for (uInt i=0; i<newShape.size(); ++i) {
758 0 : if (i == (uInt)_fitAxis) {
759 0 : newShape[i] = 1;
760 : }
761 : else {
762 0 : newShape[i] = oldShape[oldIndex];
763 0 : ++oldIndex;
764 : }
765 : }
766 0 : fitMask.assign(fitMask.reform(newShape));
767 0 : }
768 0 : _loopOverFits(
769 : fitData, showProgress, pProgressMeter, checkMinPts,
770 : fitMask, abcissaType, fitterShape, sliceShape,
771 : myGoodPlanes
772 : );
773 0 : }
774 :
775 0 : void ImageProfileFitter::_loopOverFits(
776 : SPCIIF fitData, Bool showProgress,
777 : std::shared_ptr<ProgressMeter> progressMeter, Bool checkMinPts,
778 : const Array<Bool>& fitMask, ImageFit1D<Float>::AbcissaType abcissaType,
779 : const IPosition& fitterShape, const IPosition& sliceShape,
780 : const std::set<uInt> goodPlanes
781 : ) {
782 0 : *_getLog() << LogOrigin(_class, __func__);
783 0 : Lattice<Bool>* pFitMask = _modelImage && _modelImage->hasPixelMask()
784 0 : && _modelImage->pixelMask().isWritable()
785 0 : ? &(_modelImage->pixelMask())
786 0 : : 0;
787 0 : Lattice<Bool>* pResidMask = _residImage && _residImage->hasPixelMask()
788 0 : && _residImage->pixelMask().isWritable()
789 0 : ? &(_residImage->pixelMask())
790 0 : : 0;
791 0 : vector<IPosition> goodPos(0);
792 0 : SpectralList newEstimates = _nonPolyEstimates;
793 : ImageFit1D<Float> fitter = _sigma
794 0 : ? ImageFit1D<Float>(fitData, _sigma, _fitAxis)
795 0 : : ImageFit1D<Float>(fitData, _fitAxis);
796 0 : Bool isSpectral = _fitAxis == _subImage->coordinates().spectralAxisNumber();
797 :
798 : // calculate the abscissa values only once if they will not change
799 : // with position
800 0 : Double *divisorPtr = 0;
801 0 : Vector<Double> abscissaValues(0);
802 0 : Bool fitSuccess = False;
803 0 : if (isSpectral) {
804 0 : abscissaValues = fitter.makeAbscissa(abcissaType, True, 0);
805 0 : if (_isSpectralIndex) {
806 0 : _setAbscissaDivisorIfNecessary(abscissaValues);
807 0 : if (_abscissaDivisor != 1) {
808 0 : divisorPtr = &_abscissaDivisor;
809 0 : abscissaValues /= _abscissaDivisor;
810 0 : if (_nLTPCoeffs > 0) {
811 0 : abscissaValues = log(abscissaValues);
812 : }
813 : }
814 : }
815 : }
816 0 : std::unique_ptr<const PolynomialSpectralElement> polyEl;
817 0 : if (_polyOrder >= 0) {
818 0 : polyEl.reset(new PolynomialSpectralElement(Vector<Double>(_polyOrder + 1, 0)));
819 0 : if (newEstimates.nelements() > 0) {
820 0 : newEstimates.add(*polyEl);
821 : }
822 : }
823 0 : uInt nOrigComps = newEstimates.nelements();
824 0 : Array<Double> (*xfunc)(const Array<Double>&) = 0;
825 0 : Array<Double> (*yfunc)(const Array<Double>&) = 0;
826 0 : Bool abscissaSet = ! abscissaValues.empty();
827 0 : if (_nLTPCoeffs > 0) {
828 0 : if (! abscissaSet) {
829 0 : xfunc = casacore::log;
830 : }
831 0 : yfunc = casacore::log;
832 : }
833 0 : if (abscissaSet) {
834 0 : fitter.setAbscissa(abscissaValues);
835 : //abscissaSet = False;
836 : }
837 0 : IPosition inTileShape = fitData->niceCursorShape();
838 0 : TiledLineStepper stepper (fitData->shape(), inTileShape, _fitAxis);
839 0 : RO_MaskedLatticeIterator<Float> inIter(*fitData, stepper);
840 0 : uInt nProfiles = 0;
841 0 : Bool hasXMask = ! goodPlanes.empty();
842 0 : Bool hasNonPolyEstimates = _nonPolyEstimates.nelements() > 0;
843 0 : Bool updateOutput = _modelImage || _residImage;
844 0 : Bool storeGoodPos = hasNonPolyEstimates && ! _fitters.empty();
845 0 : for (inIter.reset(); ! inIter.atEnd(); ++inIter, ++nProfiles) {
846 0 : if (showProgress && /*nProfiles % mark == 0 &&*/ nProfiles > 0) {
847 0 : progressMeter->update(Double(nProfiles));
848 : }
849 0 : const IPosition& curPos = inIter.position();
850 0 : if (checkMinPts && ! fitMask(curPos)) {
851 0 : continue;
852 : }
853 0 : ++_nAttempted;
854 0 : fitter.clearList();
855 0 : if (abscissaSet) {
856 0 : fitter.setData(curPos, yfunc);
857 : }
858 : else {
859 0 : fitter.setData(
860 : curPos, abcissaType, True, divisorPtr, xfunc, yfunc
861 : );
862 : }
863 0 : Bool canFit = _setFitterElements(
864 : fitter, newEstimates, polyEl, goodPos,
865 : fitterShape, curPos, nOrigComps
866 : );
867 0 : if (canFit) {
868 0 : if (hasXMask) {
869 0 : fitter.setXMask(goodPlanes, True);
870 : }
871 : try {
872 0 : fitSuccess = fitter.fit();
873 0 : if (fitSuccess) {
874 0 : if (fitter.converged()) {
875 0 : _flagFitterIfNecessary(fitter);
876 0 : ++_nConverged;
877 : }
878 0 : fitSuccess = fitter.isValid();
879 0 : if (fitSuccess) {
880 0 : ++_nValid;
881 0 : if (storeGoodPos) {
882 0 : goodPos.push_back(curPos);
883 : }
884 : }
885 : }
886 : }
887 0 : catch (const AipsError& x) {
888 0 : fitSuccess = False;
889 0 : }
890 : }
891 : else {
892 0 : fitSuccess = False;
893 : }
894 0 : if (fitter.succeeded()) {
895 0 : ++_nSucceeded;
896 : }
897 0 : if (_storeFits) {
898 0 : _fitters(curPos).reset(new ProfileFitResults(fitter));
899 : }
900 0 : if (updateOutput) {
901 0 : _updateModelAndResidual(
902 : fitSuccess, fitter, sliceShape,
903 : curPos, pFitMask, pResidMask
904 : );
905 : }
906 0 : }
907 0 : }
908 :
909 0 : void ImageProfileFitter::_updateModelAndResidual(
910 : Bool fitOK, const ImageFit1D<Float>& fitter,
911 : const IPosition& sliceShape, const IPosition& curPos,
912 : Lattice<Bool>* const &pFitMask,
913 : Lattice<Bool>* const &pResidMask
914 : ) const {
915 0 : static const Array<Float> failData(sliceShape, NAN);
916 0 : static const Array<Bool> failMask(sliceShape, False);
917 : Array<Bool> resultMask = fitOK
918 0 : ? fitter.getDataMask().reform(sliceShape)
919 0 : : failMask;
920 0 : if (_modelImage) {
921 0 : _modelImage->putSlice (
922 0 : (fitOK ? fitter.getFit().reform(sliceShape) : failData),
923 : curPos
924 : );
925 0 : if (pFitMask) {
926 0 : pFitMask->putSlice(resultMask, curPos);
927 : }
928 : }
929 0 : if (_residImage) {
930 0 : _residImage->putSlice (
931 0 : (fitOK ? fitter.getResidual().reform(sliceShape) : failData),
932 : curPos
933 : );
934 0 : if (pResidMask) {
935 0 : pResidMask->putSlice(resultMask, curPos);
936 : }
937 : }
938 0 : }
939 :
940 0 : Bool ImageProfileFitter::_setFitterElements(
941 : ImageFit1D<Float>& fitter, SpectralList& newEstimates,
942 : const std::unique_ptr<const PolynomialSpectralElement>& polyEl,
943 : const std::vector<IPosition>& goodPos,
944 : const IPosition& fitterShape, const IPosition& curPos,
945 : uInt nOrigComps
946 : ) const {
947 0 : if (_nonPolyEstimates.nelements() == 0) {
948 0 : if (_nGaussSinglets > 0) {
949 0 : fitter.setGaussianElements (_nGaussSinglets);
950 0 : uInt ng = fitter.getList(False).nelements();
951 0 : if (ng != _nGaussSinglets && ! _haveWarnedAboutGuessingGaussians) {
952 0 : *this->_getLog() << LogOrigin(getClass(), __func__) << LogIO::WARN;
953 0 : if (ng == 0) {
954 0 : *this->_getLog() << "Unable to estimate "
955 0 : << "parameters for any Gaussian singlets. ";
956 : }
957 : else {
958 0 : *this->_getLog() << "Only able to estimate parameters for " << ng
959 0 : << " Gaussian singlets. ";
960 : }
961 0 : *this->_getLog() << "If you really want "
962 0 : << _nGaussSinglets << " Gaussian singlets to be fit, "
963 0 : << "you should specify initial parameter estimates for all of them";
964 0 : if (_multiFit) {
965 0 : *this->_getLog() << " (additional warnings of this type during "
966 0 : "this run will not be logged)";
967 : }
968 0 : *this->_getLog() << "." << LogIO::POST;
969 0 : _haveWarnedAboutGuessingGaussians = True;
970 : }
971 : }
972 0 : if (polyEl.get()) {
973 0 : fitter.addElement(*polyEl);
974 : }
975 : else {
976 0 : if (fitter.getList(False).nelements() == 0) {
977 0 : return False;
978 : }
979 : }
980 : }
981 : else {
982 : // user supplied initial estimates
983 0 : if (goodPos.size() > 0) {
984 0 : IPosition nearest;
985 0 : Int minDist2 = 0;
986 0 : for (
987 0 : IPosition::const_iterator iter=fitterShape.begin();
988 0 : iter!=fitterShape.end(); iter++
989 : ) {
990 0 : minDist2 += *iter * *iter;
991 : }
992 0 : for (
993 0 : vector<IPosition>::const_reverse_iterator iter=goodPos.rbegin();
994 0 : iter != goodPos.rend(); iter++
995 : ) {
996 0 : IPosition diff = curPos - *iter;
997 0 : Int dist2 = 0;
998 0 : Bool larger = False;
999 0 : for (
1000 0 : IPosition::const_iterator ipositer=diff.begin();
1001 0 : ipositer!=diff.end(); ipositer++
1002 : ) {
1003 0 : dist2 += *ipositer * *ipositer;
1004 0 : if(dist2 >= minDist2) {
1005 0 : larger = True;
1006 0 : break;
1007 : }
1008 : }
1009 0 : if (
1010 0 : _fitters(*iter)->getList().nelements() == nOrigComps
1011 0 : && ! larger
1012 : ) {
1013 0 : minDist2 = dist2;
1014 0 : nearest = *iter;
1015 0 : if (minDist2 == 1) {
1016 : // can't get any nearer than this
1017 0 : break;
1018 : }
1019 : }
1020 0 : }
1021 0 : newEstimates = _fitters(nearest)->getList();
1022 0 : }
1023 0 : fitter.setElements(newEstimates);
1024 : }
1025 0 : return True;
1026 : }
1027 :
1028 0 : void ImageProfileFitter::_setAbscissaDivisorIfNecessary(
1029 : const Vector<Double>& abscissaValues
1030 : ) {
1031 0 : if (_abscissaDivisor == 0) {
1032 0 : setAbscissaDivisor(1);
1033 0 : if (abscissaValues.size() > 0) {
1034 0 : Double minAbs = min(abs(abscissaValues));
1035 0 : Double maxAbs = max(abs(abscissaValues));
1036 0 : Double l = (Int)log10(sqrt(minAbs*maxAbs));
1037 0 : Double p = std::pow(10.0, l);
1038 0 : setAbscissaDivisor(p);
1039 : }
1040 : }
1041 0 : if (_abscissaDivisor != 1) {
1042 0 : *_getLog() << LogIO::NORMAL << "Dividing abscissa values by "
1043 0 : << _abscissaDivisor << " before fitting" << LogIO::POST;
1044 : }
1045 0 : }
1046 :
1047 0 : void ImageProfileFitter::_flagFitterIfNecessary(
1048 : ImageFit1D<Float>& fitter
1049 : ) const {
1050 0 : Bool checkComps = _goodAmpRange || _goodCenterRange
1051 0 : || _goodFWHMRange;
1052 0 : SpectralList solutions = fitter.getList(True);
1053 0 : for (uInt i=0; i<solutions.nelements(); ++i) {
1054 0 : if (
1055 0 : anyTrue(isNaN(solutions[i]->get()))
1056 0 : || anyTrue(isNaN(solutions[i]->getError()))
1057 : ) {
1058 0 : fitter.invalidate();
1059 0 : return;
1060 : }
1061 0 : if (checkComps) {
1062 0 : switch (solutions[i]->getType()) {
1063 0 : case SpectralElement::GAUSSIAN:
1064 : // allow fall through
1065 : case SpectralElement::LORENTZIAN: {
1066 0 : if (
1067 0 : ! _isPCFSolutionOK(
1068 0 : dynamic_cast<
1069 : const PCFSpectralElement*
1070 0 : >(solutions[i])
1071 : )
1072 : ) {
1073 0 : fitter.invalidate();
1074 0 : return;
1075 : }
1076 0 : break;
1077 : }
1078 0 : case SpectralElement::GMULTIPLET: {
1079 0 : const GaussianMultipletSpectralElement *gm = dynamic_cast<
1080 : const GaussianMultipletSpectralElement*
1081 0 : >(solutions[i]);
1082 0 : Vector<GaussianSpectralElement> gse(gm->getGaussians());
1083 0 : for (uInt j=0; j<gse.size(); j++) {
1084 0 : if (! _isPCFSolutionOK(&gse[i])) {
1085 0 : fitter.invalidate();
1086 0 : return;
1087 : }
1088 : }
1089 0 : break;
1090 0 : }
1091 0 : default:
1092 0 : continue;
1093 0 : }
1094 : }
1095 : }
1096 0 : }
1097 :
1098 0 : Bool ImageProfileFitter::_isPCFSolutionOK(
1099 : const PCFSpectralElement *const &pcf
1100 : ) const {
1101 0 : if (_goodAmpRange) {
1102 0 : Double amp = pcf->getAmpl();
1103 0 : if (
1104 0 : amp < _goodAmpRange->first
1105 0 : || amp > _goodAmpRange->second
1106 0 : || fabs(pcf->getAmplErr()/amp) > 100
1107 : ) {
1108 0 : return False;
1109 : }
1110 : }
1111 0 : if (_goodCenterRange) {
1112 0 : Double center = pcf->getCenter();
1113 0 : if (
1114 0 : center < _goodCenterRange->first
1115 0 : || center > _goodCenterRange->second
1116 : ) {
1117 0 : return False;
1118 : }
1119 : }
1120 0 : if (_goodFWHMRange) {
1121 0 : Double fwhm = pcf->getFWHM();
1122 0 : if (
1123 0 : fwhm < _goodFWHMRange->first
1124 0 : || fwhm > _goodFWHMRange->second
1125 0 : || fabs(pcf->getFWHMErr()/fwhm) > 100
1126 : ) {
1127 0 : return False;
1128 : }
1129 : }
1130 0 : return True;
1131 : }
1132 :
1133 0 : const Array<std::shared_ptr<ProfileFitResults> >& ImageProfileFitter::getFitters() const{
1134 0 : return _fitters;
1135 : }
1136 :
1137 : }
|
