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 70 : 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 70 : ) : ImageTask<Float>(
60 : image, region, regionPtr, box, chans, stokes,
61 : mask, "", False
62 : ),
63 70 : _residual(), _model(), _xUnit(), _centerName(),
64 70 : _centerErrName(), _fwhmName(), _fwhmErrName(),
65 70 : _ampName(), _ampErrName(), _integralName(),
66 140 : _integralErrName(), _plpName(), _plpErrName(), _sigmaName(),
67 70 : _abscissaDivisorForDisplay("1"), _multiFit(False),
68 70 : /*_deleteImageOnDestruct(False),*/ _logResults(True),
69 70 : _isSpectralIndex(False), _createResid(False), _overwrite(overwrite),
70 70 : _storeFits(True),
71 70 : _polyOrder(-1), _fitAxis(axis), _nGaussSinglets(ngauss),
72 70 : _nGaussMultiplets(0), _nLorentzSinglets(0), _nPLPCoeffs(0),
73 70 : _nLTPCoeffs(0), _minGoodPoints(1), _nProfiles(0), _nAttempted(0), _nSucceeded(0),
74 70 : _nConverged(0), _nValid(0), _results(Record()), _nonPolyEstimates(),
75 70 : _goodAmpRange(), _goodCenterRange(), _goodFWHMRange(),
76 280 : _sigma(), _abscissaDivisor(1.0), _residImage(), _goodPlanes() {
77 70 : *_getLog() << LogOrigin(_class, __func__);
78 70 : _construct();
79 64 : _finishConstruction();
80 273 : }
81 :
82 5 : 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 5 : ) : ImageTask<Float>(
89 : image, region, regionPtr, box, chans, stokes,
90 : mask, "", False
91 : ),
92 5 : _residual(), _model(), _xUnit(), _centerName(),
93 5 : _centerErrName(), _fwhmName(), _fwhmErrName(),
94 5 : _ampName(), _ampErrName(), _integralName(),
95 10 : _integralErrName(), _plpName(), _plpErrName(), _sigmaName(),
96 5 : _abscissaDivisorForDisplay("1"), _multiFit(False),
97 5 : /*_deleteImageOnDestruct(False),*/ _logResults(True),
98 5 : _isSpectralIndex(False), _createResid(False), _overwrite(overwrite),
99 5 : _storeFits(True),
100 5 : _polyOrder(-1), _fitAxis(axis), _nGaussSinglets(0),
101 5 : _nGaussMultiplets(0), _nLorentzSinglets(0), _nPLPCoeffs(0),
102 5 : _nLTPCoeffs(0), _minGoodPoints(1), _nProfiles(0), _nAttempted(0), _nSucceeded(0),
103 5 : _nConverged(0), _nValid(0), _results(Record()), _nonPolyEstimates(),
104 5 : _goodAmpRange(), _goodCenterRange(), _goodFWHMRange(),
105 20 : _sigma(), _abscissaDivisor(1.0), _residImage(), _goodPlanes() {
106 5 : *_getLog() << LogOrigin(_class, __func__);
107 5 : ThrowIf(estimatesFilename.empty(), "Estimates filename cannot be empty");
108 5 : ProfileFitterEstimatesFileParser parser(estimatesFilename);
109 5 : _nonPolyEstimates = parser.getEstimates();
110 5 : _nGaussSinglets = _nonPolyEstimates.nelements();
111 10 : *_getLog() << LogIO::NORMAL << "Number of gaussian singlets to fit found to be "
112 : <<_nGaussSinglets << " in estimates file " << estimatesFilename
113 5 : << LogIO::POST;
114 5 : _construct();
115 5 : _finishConstruction();
116 5 : }
117 :
118 21 : 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 21 : ) : ImageTask<Float>(
125 : image, region, regionPtr, box, chans, stokes,
126 : mask, "", False
127 : ),
128 21 : _residual(), _model(), _xUnit(), _centerName(),
129 21 : _centerErrName(), _fwhmName(), _fwhmErrName(),
130 21 : _ampName(), _ampErrName(), _integralName(),
131 42 : _integralErrName(), _plpName(), _plpErrName(), _sigmaName(),
132 21 : _abscissaDivisorForDisplay("1"), _multiFit(False),
133 21 : /* _deleteImageOnDestruct(False),*/ _logResults(True),
134 21 : _isSpectralIndex(False), _createResid(False), _overwrite(overwrite),
135 21 : _storeFits(True),
136 21 : _polyOrder(-1), _fitAxis(axis), _nGaussSinglets(0),
137 21 : _nGaussMultiplets(0), _nLorentzSinglets(0), _nPLPCoeffs(0),
138 21 : _nLTPCoeffs(0), _minGoodPoints(1), _nProfiles(0), _nAttempted(0), _nSucceeded(0),
139 21 : _nConverged(0), _nValid(0), _results(Record()), _nonPolyEstimates(),
140 21 : _goodAmpRange(), _goodCenterRange(), _goodFWHMRange(),
141 84 : _sigma(), _abscissaDivisor(1.0), _residImage(), _goodPlanes() {
142 21 : *_getLog() << LogOrigin(_class, __func__);
143 21 : ThrowIf(
144 : spectralList.nelements() == 0,
145 : "spectralList cannot be empty"
146 : );
147 21 : _nonPolyEstimates = spectralList;
148 21 : _nGaussSinglets = 0;
149 21 : _nGaussMultiplets = 0;
150 45 : for (uInt i=0; i<_nonPolyEstimates.nelements(); i++) {
151 24 : SpectralElement::Types myType = _nonPolyEstimates[i]->getType();
152 24 : switch(myType) {
153 4 : case SpectralElement::GAUSSIAN:
154 4 : _nGaussSinglets++;
155 4 : break;
156 2 : case SpectralElement::GMULTIPLET:
157 2 : _nGaussMultiplets++;
158 2 : break;
159 2 : case SpectralElement::LORENTZIAN:
160 2 : _nLorentzSinglets++;
161 2 : break;
162 11 : case SpectralElement::POWERLOGPOLY:
163 11 : 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 11 : _nPLPCoeffs = _nonPolyEstimates[i]->get().size();
169 11 : break;
170 5 : case SpectralElement::LOGTRANSPOLY:
171 5 : 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 5 : _nLTPCoeffs = _nonPolyEstimates[i]->get().size();
177 5 : 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 24 : if (_nPLPCoeffs > 0) {
188 22 : *_getLog() << LogIO::NORMAL << "Will fit a single power logarithmic polynomial "
189 11 : << " from provided spectral element list" << LogIO::POST;
190 : }
191 13 : else if (_nLTPCoeffs > 0) {
192 10 : *_getLog() << LogIO::NORMAL << "Will fit a single logarithmic transformed polynomial "
193 5 : << " from provided spectral element list" << LogIO::POST;
194 : }
195 : else {
196 8 : if (_nGaussSinglets > 0) {
197 8 : *_getLog() << LogIO::NORMAL << "Number of Gaussian singlets to fit found to be "
198 : << _nGaussSinglets << " from provided spectral element list"
199 4 : << LogIO::POST;
200 : }
201 8 : if (_nGaussMultiplets > 0) {
202 4 : *_getLog() << LogIO::NORMAL << "Number of Gaussian multiplets to fit found to be "
203 : << _nGaussMultiplets << " from provided spectral element list"
204 2 : << LogIO::POST;
205 : }
206 8 : if (_nLorentzSinglets > 0) {
207 4 : *_getLog() << LogIO::NORMAL << "Number of lorentzian singlets to fit found to be "
208 : << _nLorentzSinglets << " from provided spectral element list"
209 2 : << LogIO::POST;
210 : }
211 : }
212 : }
213 21 : _construct();
214 21 : _finishConstruction();
215 21 : }
216 :
217 175 : ImageProfileFitter::~ImageProfileFitter() {}
218 :
219 89 : 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 89 : _checkNGaussAndPolyOrder();
223 88 : LogOrigin logOrigin(_class, __func__);
224 88 : *_getLog() << logOrigin;
225 88 : std::unique_ptr<ImageInterface<Float> > originalSigma;
226 : {
227 264 : _subImage = SubImageFactory<Float>::createSubImageRO(
228 177 : *_getImage(), *_getRegion(), _getMask(), 0,
229 177 : AxesSpecifier(), _getStretch()
230 87 : );
231 87 : uInt nUnknowns = _nUnknowns();
232 94 : 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 86 : if (_sigma.get()) {
240 48 : if (! _sigmaName.empty()) {
241 48 : originalSigma.reset(_sigma->cloneII());
242 : }
243 : std::shared_ptr<const SubImage<Float> > sigmaSubImage = SubImageFactory<Float>::createSubImageRO(
244 96 : *_sigma, *_getRegion(), _getMask(), 0, AxesSpecifier(), _getStretch()
245 96 : );
246 96 : _sigma.reset(
247 : new TempImage<Float>(
248 144 : sigmaSubImage->shape(), sigmaSubImage->coordinates()
249 48 : )
250 : );
251 48 : _sigma->put(sigmaSubImage->get());
252 48 : }
253 : }
254 86 : *_getLog() << logOrigin;
255 3 : _storeFits = doDetailedResults || ! _centerName.empty()
256 3 : || ! _centerErrName.empty() || ! _fwhmName.empty()
257 3 : || ! _fwhmErrName.empty() || ! _ampName.empty()
258 3 : || ! _ampErrName.empty() || ! _integralName.empty()
259 3 : || ! _integralErrName.empty()
260 3 : || ! _plpName.empty() || ! _plpErrName.empty()
261 89 : || ! _ltpName.empty() || ! _ltpErrName.empty();
262 : try {
263 86 : if (! _multiFit) {
264 : ImageCollapser<Float> collapser(
265 76 : _subImage, IPosition(1, _fitAxis), True,
266 : ImageCollapserData::MEAN, "", True
267 114 : );
268 38 : SPIIF x = collapser.collapse();
269 : // _subImage needs to be a SubImage<Float> object
270 76 : _subImage = SubImageFactory<Float>::createSubImageRO(
271 76 : *x, Record(), "", _getLog().get(),
272 76 : AxesSpecifier(), False
273 38 : );
274 38 : if (_sigma.get()) {
275 48 : Array<Bool> sigmaMask = _sigma->get() != Array<Float>(_sigma->shape(), 0.0f);
276 24 : if (anyTrue(! sigmaMask)) {
277 6 : if (_sigma->hasPixelMask()) {
278 0 : sigmaMask = sigmaMask && _sigma->pixelMask().get();
279 : }
280 : else {
281 6 : _sigma->makeMask("sigmamask", True, True, False);
282 : }
283 6 : _sigma->pixelMask().put(sigmaMask);
284 : }
285 : ImageCollapser<Float> collapsedSigma(
286 48 : _sigma, IPosition(1, _fitAxis), True,
287 : ImageCollapserData::MEAN, "", True
288 72 : );
289 24 : SPIIF collapsed = collapsedSigma.collapse();
290 24 : std::shared_ptr<TempImage<Float> >ctmp = std::dynamic_pointer_cast<TempImage<Float> >(collapsed);
291 24 : ThrowIf(! ctmp, "Dynamic cast failed");
292 24 : _sigma = ctmp;
293 24 : }
294 38 : }
295 86 : _fitallprofiles();
296 86 : *_getLog() << logOrigin;
297 : }
298 0 : catch (const AipsError& x) {
299 0 : ThrowCc("Exception during fit: " + x.getMesg());
300 0 : }
301 : ImageProfileFitterResults resultHandler(
302 86 : _getLog(), _getImage()->coordinates(), &_fitters,
303 86 : _nonPolyEstimates, _subImage, _polyOrder,
304 : _fitAxis, _nGaussSinglets, _nGaussMultiplets,
305 86 : _nLorentzSinglets, _nPLPCoeffs, _nLTPCoeffs, _logResults,
306 172 : _multiFit, _getLogFile(), _xUnit, _summaryHeader()
307 258 : );
308 86 : addHistory(
309 : logOrigin,
310 172 : resultHandler.logSummary(
311 : _nProfiles, _nAttempted, _nSucceeded, _nConverged, _nValid
312 : )
313 : );
314 86 : if (_nPLPCoeffs > 0) {
315 11 : resultHandler.setPLPName(_plpName);
316 11 : resultHandler.setPLPErrName(_plpErrName);
317 11 : resultHandler.setPLPDivisor(_abscissaDivisorForDisplay);
318 : }
319 75 : else if (_nLTPCoeffs > 0) {
320 5 : resultHandler.setLTPName(_ltpName);
321 5 : resultHandler.setLTPErrName(_ltpErrName);
322 5 : resultHandler.setPLPDivisor(_abscissaDivisorForDisplay);
323 : }
324 70 : else if (_nGaussSinglets + _nGaussMultiplets + _nLorentzSinglets > 0) {
325 67 : resultHandler.setAmpName(_ampName);
326 67 : resultHandler.setAmpErrName(_ampErrName);
327 67 : resultHandler.setCenterName(_centerName);
328 67 : resultHandler.setCenterErrName(_centerErrName);
329 67 : resultHandler.setFWHMName(_fwhmName);
330 67 : resultHandler.setFWHMErrName(_fwhmErrName);
331 67 : resultHandler.setIntegralName(_integralName);
332 67 : resultHandler.setIntegralErrName(_integralErrName);
333 : }
334 86 : if (doDetailedResults) {
335 83 : resultHandler.createResults();
336 83 : _results = resultHandler.getResults();
337 : }
338 86 : resultHandler.writeImages(_nConverged > 0);
339 86 : if (_modelImage) {
340 5 : _modelImage = _prepareOutputImage(*_modelImage, _model, _overwrite, True);
341 : }
342 86 : if (_residImage) {
343 5 : _residImage = _prepareOutputImage(*_residImage, _residual, _overwrite, True);
344 : }
345 86 : if (originalSigma && ! _sigmaName.empty()) {
346 48 : _prepareOutputImage(*originalSigma, _sigmaName, True, True);
347 : }
348 172 : return _results;
349 90 : }
350 :
351 87 : uInt ImageProfileFitter::_nUnknowns() const {
352 87 : uInt n = 0;
353 87 : if (_polyOrder >= 0) {
354 8 : n += _polyOrder + 1;
355 : }
356 87 : if (_nGaussSinglets > 0) {
357 65 : n += 3*_nGaussSinglets;
358 : }
359 87 : uInt nel = _nonPolyEstimates.nelements();
360 87 : if (n == 0) {
361 19 : return n;
362 : }
363 79 : for (uInt i=0; i<nel; ++i) {
364 11 : const SpectralElement *const x = _nonPolyEstimates[i];
365 11 : Vector<Bool> fixed = x->fixed();
366 11 : Vector<Bool>::const_iterator iter = fixed.begin();
367 11 : Vector<Bool>::const_iterator end = fixed.end();
368 44 : while (iter != end) {
369 33 : if (*iter) {
370 1 : --n;
371 1 : if (n == 0) {
372 0 : return n;
373 : }
374 : }
375 33 : ++iter;
376 : }
377 11 : }
378 68 : return n;
379 : }
380 :
381 :
382 8 : void ImageProfileFitter::setPolyOrder(Int p) {
383 8 : ThrowIf(p < 0,"A polynomial cannot have a negative order");
384 8 : ThrowIf(
385 : _nPLPCoeffs > 0,
386 : "Cannot simultaneously fit a polynomial and "
387 : "a power logarithmic polynomial."
388 : );
389 8 : ThrowIf(
390 : _nLTPCoeffs > 0,
391 : "Cannot simultaneously fit a polynomial and "
392 : "a logarithmic transformed polynomial"
393 : );
394 8 : _polyOrder = p;
395 8 : }
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 24 : void ImageProfileFitter::setSigma(const Array<Float>& sigma) {
416 24 : std::unique_ptr<TempImage<Float> > temp;
417 24 : if (sigma.ndim() == _getImage()->ndim()) {
418 48 : temp.reset(new TempImage<Float>(
419 32 : sigma.shape(), _getImage()->coordinates())
420 : );
421 : }
422 8 : else if (sigma.ndim() == 1) {
423 8 : SpectralCoordinate sp;
424 8 : CoordinateSystem csys;
425 8 : csys.addCoordinate(sp);
426 8 : temp.reset(new TempImage<Float>(sigma.shape(), csys));
427 8 : }
428 24 : temp->put(sigma);
429 24 : setSigma(temp.get());
430 24 : }
431 :
432 48 : void ImageProfileFitter::setSigma(const ImageInterface<Float>* const &sigma) {
433 48 : if (anyTrue(sigma->get() < Array<Float>(sigma->shape(), 0.0f))) {
434 0 : *_getLog() << "All sigma values must be non-negative" << LogIO::EXCEPTION;
435 : }
436 48 : Float mymax = fabs(max(sigma->get()));
437 48 : if (
438 144 : sigma->ndim() == _getImage()->ndim()
439 144 : && sigma->shape() == _getImage()->shape()
440 : ) {
441 16 : SPIIF clone(sigma->cloneII());
442 16 : _sigma = std::dynamic_pointer_cast<TempImage<Float> >(clone);
443 16 : if (! _sigma) {
444 : SPIIF x = SubImageFactory<Float>::createImage(
445 16 : *sigma, "", Record(), "", False, False ,True, False
446 24 : );
447 8 : if (x) {
448 8 : _sigma = std::dynamic_pointer_cast<TempImage<Float> >(x);
449 8 : ThrowIf(
450 : ! _sigma,
451 : "Unable to create temporary weights image"
452 : );
453 : }
454 8 : }
455 16 : if (mymax != 1) {
456 4 : _sigma->put(_sigma->get()/mymax);
457 : }
458 16 : }
459 32 : else if (
460 96 : sigma->ndim() == _getImage()->ndim()
461 96 : || sigma->ndim() == 1
462 : ) {
463 32 : if (sigma->ndim() == _getImage()->ndim()) {
464 16 : IPosition expShape(_getImage()->ndim(), 1);
465 16 : expShape[_fitAxis] = _getImage()->shape()[_fitAxis];
466 16 : 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 16 : }
473 16 : else if (sigma->ndim() == 1) {
474 16 : 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 32 : IPosition dataToInsertShape(_getImage()->ndim(), 1);
481 32 : dataToInsertShape[_fitAxis] = _getImage()->shape()[_fitAxis];
482 32 : _sigma.reset(new TempImage<Float>(_getImage()->shape(), _getImage()->coordinates()));
483 64 : Array<Float> dataToInsert(IPosition(1, _getImage()->shape()[_fitAxis]));
484 32 : dataToInsert = sigma->get(sigma->ndim() == _getImage()->ndim());
485 : // normalize
486 32 : if (mymax != 1) {
487 8 : dataToInsert /= mymax;
488 : }
489 32 : Array<Float> sigmaData = _sigma->get();
490 32 : ArrayIterator<Float> iter(sigmaData, IPosition(1, _fitAxis), True);
491 2624 : while(! iter.pastEnd()) {
492 2592 : iter.array() = dataToInsert;
493 2592 : iter.next();
494 : }
495 32 : _sigma->put(sigmaData);
496 32 : }
497 : else {
498 0 : ThrowCc("Illegal shape of standard deviation image");
499 : }
500 48 : if (! _sigma->coordinates().near(_getImage()->coordinates())) {
501 8 : _sigma->setCoordinateInfo(_getImage()->coordinates());
502 : }
503 48 : }
504 :
505 0 : Record ImageProfileFitter::getResults() const {
506 0 : return _results;
507 : }
508 :
509 45 : void ImageProfileFitter::setAbscissaDivisor(Double d) {
510 45 : 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 45 : _abscissaDivisor = d;
517 0 : _abscissaDivisorForDisplay = String::toString(d)
518 45 : + _getImage()->coordinates().worldAxisUnits()[_fitAxis];
519 45 : }
520 :
521 1 : void ImageProfileFitter::setAbscissaDivisor(const Quantity& q) {
522 2 : String fitAxisUnit = _getImage()->coordinates().worldAxisUnits()[_fitAxis];
523 1 : ThrowIf(
524 : ! q.isConform(fitAxisUnit),
525 : "Abscissa divisor unit " + q.getUnit()
526 : + " is not consistent with fit axis unit"
527 : );
528 1 : 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 1 : _abscissaDivisor = q.getValue(fitAxisUnit);
535 1 : _abscissaDivisorForDisplay = String::toString(q);
536 1 : }
537 :
538 96 : std::vector<OutputDestinationChecker::OutputStruct> ImageProfileFitter::_getOutputStruct() {
539 96 : std::vector<OutputDestinationChecker::OutputStruct> outputs;
540 96 : 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 96 : 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 96 : return outputs;
557 0 : }
558 :
559 89 : void ImageProfileFitter::_checkNGaussAndPolyOrder() const {
560 89 : 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 88 : }
570 :
571 90 : void ImageProfileFitter::_finishConstruction() {
572 90 : LogOrigin logOrigin(_class, __func__);
573 90 : _isSpectralIndex = _nPLPCoeffs + _nLTPCoeffs > 0;
574 97 : 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 89 : if (_fitAxis < 0) {
581 20 : 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 20 : _fitAxis = _getImage()->coordinates().spectralAxisNumber();
588 40 : *_getLog() << LogIO::NORMAL << "Using spectral axis (axis " << _fitAxis
589 20 : << ") as fit axis" << LogIO::POST;
590 : }
591 : }
592 : //this->_setSupportsLogfile(true);
593 90 : }
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 86 : void ImageProfileFitter::_fitallprofiles() {
645 86 : *_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 86 : if (
650 172 : ! _model.empty()
651 91 : && ! (
652 182 : _modelImage = SubImageFactory<Float>::createImage(
653 96 : *_subImage, "", Record(), "",
654 5 : False, _overwrite ,False, False, True
655 5 : )
656 5 : )
657 : ) {
658 0 : *_getLog() << LogIO::WARN << "Failed to create model image" << LogIO::POST;
659 : }
660 86 : if (
661 81 : (! _residual.empty() || _createResid)
662 172 : && ! (
663 182 : _residImage = SubImageFactory<Float>::createImage(
664 96 : *_subImage, "", Record(), "",
665 5 : False, _overwrite ,False, False, True
666 5 : )
667 5 : )
668 : ) {
669 0 : *_getLog() << LogIO::WARN << "Failed to create residual image" << LogIO::POST;
670 : }
671 86 : Bool showProgress = True;
672 86 : _fitProfiles(showProgress);
673 86 : }
674 :
675 86 : void ImageProfileFitter::_fitProfiles(Bool showProgress) {
676 86 : IPosition inShape = _subImage->shape();
677 86 : if (_modelImage) {
678 5 : AlwaysAssert(inShape.isEqual(_modelImage->shape()), AipsError);
679 : }
680 86 : if (_residImage) {
681 5 : AlwaysAssert(inShape.isEqual(_residImage->shape()), AipsError);
682 : }
683 86 : const auto nDim = _subImage->ndim();
684 86 : IPosition sliceShape(nDim, 1);
685 86 : sliceShape(_fitAxis) = inShape(_fitAxis);
686 86 : Array<Float> resultData(sliceShape);
687 86 : Array<Bool> resultMask(sliceShape);
688 86 : String doppler = "";
689 86 : auto csys = _subImage->coordinates();
690 86 : _xUnit = csys.spectralCoordinate().worldAxisUnits()[0];
691 86 : if (
692 70 : ! _isSpectralIndex && _fitAxis == _subImage->coordinates().spectralAxisNumber()
693 156 : && _subImage->coordinates().spectralCoordinate().restFrequency() > 0
694 : ) {
695 70 : SpectralCoordinate specCoord = csys.spectralCoordinate();
696 70 : _xUnit = specCoord.velocityUnit();
697 : doppler = MDoppler::showType(
698 : specCoord.velocityDoppler()
699 70 : );
700 70 : }
701 86 : String errMsg;
702 : ImageFit1D<Float>::AbcissaType abcissaType;
703 86 : auto abscissaUnits = _isSpectralIndex ? "native" : "pix";
704 86 : ThrowIf(
705 : ! ImageFit1D<Float>::setAbcissaState(
706 : errMsg, abcissaType, csys, abscissaUnits, doppler, _fitAxis
707 : ), errMsg
708 : );
709 86 : IPosition fitterShape = inShape;
710 86 : fitterShape[_fitAxis] = 1;
711 86 : if (_storeFits) {
712 83 : _fitters.resize(fitterShape);
713 : }
714 86 : _nProfiles = fitterShape.product();
715 86 : std::shared_ptr<ProgressMeter> pProgressMeter;
716 86 : if (showProgress) {
717 86 : ostringstream oss;
718 86 : oss << "Fit profiles on axis " << _fitAxis+1;
719 172 : pProgressMeter.reset(
720 258 : new ProgressMeter(0, _nProfiles, oss.str())
721 : );
722 86 : }
723 86 : SPCIIF fitData = _subImage;
724 86 : std::set<uInt> myGoodPlanes;
725 86 : if (! _goodPlanes.empty()) {
726 1 : IPosition origin(_subImage->ndim(), 0);
727 1 : Vector<Double> world(_subImage->ndim(), 0);
728 1 : csys.toWorld(world, origin);
729 1 : const CoordinateSystem imcsys = _getImage()->coordinates();
730 1 : Int imageOff = Int(imcsys.toPixel(world)[_fitAxis] + 0.5);
731 1 : AlwaysAssert(imageOff >= 0, AipsError);
732 1 : std::vector<Int> goodPlanes(_goodPlanes.begin(), _goodPlanes.end());
733 1 : 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 1 : std::vector<Int>::iterator iter = goodPlanes.begin();
741 1 : while (iter != goodPlanes.end() && *iter < 0) {
742 0 : goodPlanes.erase(iter);
743 0 : iter = goodPlanes.begin();
744 : }
745 1 : myGoodPlanes = std::set<uInt>(goodPlanes.begin(), goodPlanes.end());
746 1 : }
747 86 : Bool checkMinPts = fitData->isMasked();
748 86 : Array<Bool> fitMask;
749 86 : if (checkMinPts) {
750 : fitMask = (
751 76 : partialNTrue(fitData->getMask(False), IPosition(1, _fitAxis))
752 114 : >= (long unsigned int) _minGoodPoints
753 38 : );
754 38 : IPosition oldShape = fitMask.shape();
755 38 : IPosition newShape(fitMask.ndim() + 1);
756 38 : uInt oldIndex = 0;
757 187 : for (uInt i=0; i<newShape.size(); ++i) {
758 149 : if (i == (uInt)_fitAxis) {
759 38 : newShape[i] = 1;
760 : }
761 : else {
762 111 : newShape[i] = oldShape[oldIndex];
763 111 : ++oldIndex;
764 : }
765 : }
766 38 : fitMask.assign(fitMask.reform(newShape));
767 38 : }
768 86 : _loopOverFits(
769 : fitData, showProgress, pProgressMeter, checkMinPts,
770 : fitMask, abcissaType, fitterShape, sliceShape,
771 : myGoodPlanes
772 : );
773 86 : }
774 :
775 86 : 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 86 : *_getLog() << LogOrigin(_class, __func__);
783 91 : Lattice<Bool>* pFitMask = _modelImage && _modelImage->hasPixelMask()
784 5 : && _modelImage->pixelMask().isWritable()
785 91 : ? &(_modelImage->pixelMask())
786 86 : : 0;
787 91 : Lattice<Bool>* pResidMask = _residImage && _residImage->hasPixelMask()
788 5 : && _residImage->pixelMask().isWritable()
789 91 : ? &(_residImage->pixelMask())
790 86 : : 0;
791 86 : vector<IPosition> goodPos(0);
792 86 : SpectralList newEstimates = _nonPolyEstimates;
793 : ImageFit1D<Float> fitter = _sigma
794 48 : ? ImageFit1D<Float>(fitData, _sigma, _fitAxis)
795 172 : : ImageFit1D<Float>(fitData, _fitAxis);
796 86 : Bool isSpectral = _fitAxis == _subImage->coordinates().spectralAxisNumber();
797 :
798 : // calculate the abscissa values only once if they will not change
799 : // with position
800 86 : Double *divisorPtr = 0;
801 86 : Vector<Double> abscissaValues(0);
802 86 : Bool fitSuccess = False;
803 86 : if (isSpectral) {
804 86 : abscissaValues = fitter.makeAbscissa(abcissaType, True, 0);
805 86 : if (_isSpectralIndex) {
806 16 : _setAbscissaDivisorIfNecessary(abscissaValues);
807 16 : if (_abscissaDivisor != 1) {
808 16 : divisorPtr = &_abscissaDivisor;
809 16 : abscissaValues /= _abscissaDivisor;
810 16 : if (_nLTPCoeffs > 0) {
811 5 : abscissaValues = log(abscissaValues);
812 : }
813 : }
814 : }
815 : }
816 86 : std::unique_ptr<const PolynomialSpectralElement> polyEl;
817 86 : if (_polyOrder >= 0) {
818 8 : polyEl.reset(new PolynomialSpectralElement(Vector<Double>(_polyOrder + 1, 0)));
819 8 : if (newEstimates.nelements() > 0) {
820 4 : newEstimates.add(*polyEl);
821 : }
822 : }
823 86 : uInt nOrigComps = newEstimates.nelements();
824 86 : Array<Double> (*xfunc)(const Array<Double>&) = 0;
825 86 : Array<Double> (*yfunc)(const Array<Double>&) = 0;
826 86 : Bool abscissaSet = ! abscissaValues.empty();
827 86 : if (_nLTPCoeffs > 0) {
828 5 : if (! abscissaSet) {
829 0 : xfunc = casacore::log;
830 : }
831 5 : yfunc = casacore::log;
832 : }
833 86 : if (abscissaSet) {
834 86 : fitter.setAbscissa(abscissaValues);
835 : //abscissaSet = False;
836 : }
837 86 : IPosition inTileShape = fitData->niceCursorShape();
838 86 : TiledLineStepper stepper (fitData->shape(), inTileShape, _fitAxis);
839 86 : RO_MaskedLatticeIterator<Float> inIter(*fitData, stepper);
840 86 : uInt nProfiles = 0;
841 86 : Bool hasXMask = ! goodPlanes.empty();
842 86 : Bool hasNonPolyEstimates = _nonPolyEstimates.nelements() > 0;
843 86 : Bool updateOutput = _modelImage || _residImage;
844 86 : Bool storeGoodPos = hasNonPolyEstimates && ! _fitters.empty();
845 808845 : for (inIter.reset(); ! inIter.atEnd(); ++inIter, ++nProfiles) {
846 808759 : if (showProgress && /*nProfiles % mark == 0 &&*/ nProfiles > 0) {
847 808673 : progressMeter->update(Double(nProfiles));
848 : }
849 808759 : const IPosition& curPos = inIter.position();
850 808759 : if (checkMinPts && ! fitMask(curPos)) {
851 263075 : continue;
852 : }
853 545684 : ++_nAttempted;
854 545684 : fitter.clearList();
855 545684 : if (abscissaSet) {
856 545684 : fitter.setData(curPos, yfunc);
857 : }
858 : else {
859 0 : fitter.setData(
860 : curPos, abcissaType, True, divisorPtr, xfunc, yfunc
861 : );
862 : }
863 545684 : Bool canFit = _setFitterElements(
864 : fitter, newEstimates, polyEl, goodPos,
865 : fitterShape, curPos, nOrigComps
866 : );
867 545684 : if (canFit) {
868 545637 : if (hasXMask) {
869 4641 : fitter.setXMask(goodPlanes, True);
870 : }
871 : try {
872 545637 : fitSuccess = fitter.fit();
873 545550 : if (fitSuccess) {
874 542204 : if (fitter.converged()) {
875 542204 : _flagFitterIfNecessary(fitter);
876 542204 : ++_nConverged;
877 : }
878 542204 : fitSuccess = fitter.isValid();
879 542204 : if (fitSuccess) {
880 542204 : ++_nValid;
881 542204 : if (storeGoodPos) {
882 537 : goodPos.push_back(curPos);
883 : }
884 : }
885 : }
886 : }
887 87 : catch (const AipsError& x) {
888 87 : fitSuccess = False;
889 87 : }
890 : }
891 : else {
892 47 : fitSuccess = False;
893 : }
894 545684 : if (fitter.succeeded()) {
895 545550 : ++_nSucceeded;
896 : }
897 545684 : if (_storeFits) {
898 16755 : _fitters(curPos).reset(new ProfileFitResults(fitter));
899 : }
900 545684 : if (updateOutput) {
901 529098 : _updateModelAndResidual(
902 : fitSuccess, fitter, sliceShape,
903 : curPos, pFitMask, pResidMask
904 : );
905 : }
906 808759 : }
907 86 : }
908 :
909 529098 : 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 529098 : static const Array<Float> failData(sliceShape, NAN);
916 529098 : static const Array<Bool> failMask(sliceShape, False);
917 : Array<Bool> resultMask = fitOK
918 1058196 : ? fitter.getDataMask().reform(sliceShape)
919 1058196 : : failMask;
920 529098 : if (_modelImage) {
921 1058138 : _modelImage->putSlice (
922 1058138 : (fitOK ? fitter.getFit().reform(sliceShape) : failData),
923 : curPos
924 : );
925 529069 : if (pFitMask) {
926 529069 : pFitMask->putSlice(resultMask, curPos);
927 : }
928 : }
929 529098 : if (_residImage) {
930 1058194 : _residImage->putSlice (
931 1058194 : (fitOK ? fitter.getResidual().reform(sliceShape) : failData),
932 : curPos
933 : );
934 529097 : if (pResidMask) {
935 529097 : pResidMask->putSlice(resultMask, curPos);
936 : }
937 : }
938 529098 : }
939 :
940 545684 : 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 545684 : if (_nonPolyEstimates.nelements() == 0) {
948 545140 : if (_nGaussSinglets > 0) {
949 16211 : fitter.setGaussianElements (_nGaussSinglets);
950 16211 : uInt ng = fitter.getList(False).nelements();
951 16211 : if (ng != _nGaussSinglets && ! _haveWarnedAboutGuessingGaussians) {
952 32 : *this->_getLog() << LogOrigin(getClass(), __func__) << LogIO::WARN;
953 32 : if (ng == 0) {
954 0 : *this->_getLog() << "Unable to estimate "
955 0 : << "parameters for any Gaussian singlets. ";
956 : }
957 : else {
958 64 : *this->_getLog() << "Only able to estimate parameters for " << ng
959 32 : << " Gaussian singlets. ";
960 : }
961 64 : *this->_getLog() << "If you really want "
962 32 : << _nGaussSinglets << " Gaussian singlets to be fit, "
963 32 : << "you should specify initial parameter estimates for all of them";
964 32 : if (_multiFit) {
965 12 : *this->_getLog() << " (additional warnings of this type during "
966 6 : "this run will not be logged)";
967 : }
968 32 : *this->_getLog() << "." << LogIO::POST;
969 32 : _haveWarnedAboutGuessingGaussians = True;
970 : }
971 : }
972 545140 : if (polyEl.get()) {
973 529010 : fitter.addElement(*polyEl);
974 : }
975 : else {
976 16130 : if (fitter.getList(False).nelements() == 0) {
977 47 : return False;
978 : }
979 : }
980 : }
981 : else {
982 : // user supplied initial estimates
983 544 : if (goodPos.size() > 0) {
984 516 : IPosition nearest;
985 516 : Int minDist2 = 0;
986 516 : for (
987 516 : IPosition::const_iterator iter=fitterShape.begin();
988 2475 : iter!=fitterShape.end(); iter++
989 : ) {
990 1959 : minDist2 += *iter * *iter;
991 : }
992 1168 : for (
993 516 : vector<IPosition>::const_reverse_iterator iter=goodPos.rbegin();
994 1684 : iter != goodPos.rend(); iter++
995 : ) {
996 1659 : IPosition diff = curPos - *iter;
997 1659 : Int dist2 = 0;
998 1659 : Bool larger = False;
999 5304 : for (
1000 1659 : IPosition::const_iterator ipositer=diff.begin();
1001 5304 : ipositer!=diff.end(); ipositer++
1002 : ) {
1003 4391 : dist2 += *ipositer * *ipositer;
1004 4391 : if(dist2 >= minDist2) {
1005 746 : larger = True;
1006 746 : break;
1007 : }
1008 : }
1009 1659 : if (
1010 3318 : _fitters(*iter)->getList().nelements() == nOrigComps
1011 3318 : && ! larger
1012 : ) {
1013 913 : minDist2 = dist2;
1014 913 : nearest = *iter;
1015 913 : if (minDist2 == 1) {
1016 : // can't get any nearer than this
1017 491 : break;
1018 : }
1019 : }
1020 1659 : }
1021 516 : newEstimates = _fitters(nearest)->getList();
1022 516 : }
1023 544 : fitter.setElements(newEstimates);
1024 : }
1025 545637 : return True;
1026 : }
1027 :
1028 16 : void ImageProfileFitter::_setAbscissaDivisorIfNecessary(
1029 : const Vector<Double>& abscissaValues
1030 : ) {
1031 16 : if (_abscissaDivisor == 0) {
1032 15 : setAbscissaDivisor(1);
1033 15 : if (abscissaValues.size() > 0) {
1034 15 : Double minAbs = min(abs(abscissaValues));
1035 15 : Double maxAbs = max(abs(abscissaValues));
1036 15 : Double l = (Int)log10(sqrt(minAbs*maxAbs));
1037 15 : Double p = std::pow(10.0, l);
1038 15 : setAbscissaDivisor(p);
1039 : }
1040 : }
1041 16 : if (_abscissaDivisor != 1) {
1042 32 : *_getLog() << LogIO::NORMAL << "Dividing abscissa values by "
1043 16 : << _abscissaDivisor << " before fitting" << LogIO::POST;
1044 : }
1045 16 : }
1046 :
1047 542204 : void ImageProfileFitter::_flagFitterIfNecessary(
1048 : ImageFit1D<Float>& fitter
1049 : ) const {
1050 1084408 : Bool checkComps = _goodAmpRange || _goodCenterRange
1051 1084408 : || _goodFWHMRange;
1052 542204 : SpectralList solutions = fitter.getList(True);
1053 1097384 : for (uInt i=0; i<solutions.nelements(); ++i) {
1054 555180 : if (
1055 1110360 : anyTrue(isNaN(solutions[i]->get()))
1056 1110360 : || anyTrue(isNaN(solutions[i]->getError()))
1057 : ) {
1058 0 : fitter.invalidate();
1059 0 : return;
1060 : }
1061 555180 : 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 542204 : }
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 : }
|
