LCOV - code coverage report
Current view: top level - imageanalysis/ImageAnalysis - ImageProfileFitter.cc (source / functions) Hit Total Coverage
Test: ctest_coverage.info Lines: 0 649 0.0 %
Date: 2023-11-02 14:27:30 Functions: 0 29 0.0 %

          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: aips2-request@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 &regionPtr,    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             :     _residual(), _model(), _xUnit(), _centerName(),
      64             :     _centerErrName(), _fwhmName(), _fwhmErrName(),
      65             :     _ampName(), _ampErrName(), _integralName(),
      66             :     _integralErrName(), _plpName(), _plpErrName(), _sigmaName(),
      67             :     _abscissaDivisorForDisplay("1"), _multiFit(False),
      68             :     /*_deleteImageOnDestruct(False),*/ _logResults(True),
      69             :     _isSpectralIndex(False), _createResid(False), _overwrite(overwrite),
      70             :     _storeFits(True),
      71             :     _polyOrder(-1), _fitAxis(axis), _nGaussSinglets(ngauss),
      72             :     _nGaussMultiplets(0), _nLorentzSinglets(0), _nPLPCoeffs(0),
      73             :     _nLTPCoeffs(0), _minGoodPoints(1), _nProfiles(0), _nAttempted(0), _nSucceeded(0),
      74             :     _nConverged(0), _nValid(0), _results(Record()), _nonPolyEstimates(),
      75             :     _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 &regionPtr,    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             :     _residual(), _model(), _xUnit(), _centerName(),
      93             :     _centerErrName(), _fwhmName(), _fwhmErrName(),
      94             :     _ampName(), _ampErrName(), _integralName(),
      95             :     _integralErrName(), _plpName(), _plpErrName(), _sigmaName(),
      96             :     _abscissaDivisorForDisplay("1"), _multiFit(False),
      97             :     /*_deleteImageOnDestruct(False),*/ _logResults(True),
      98             :     _isSpectralIndex(False), _createResid(False), _overwrite(overwrite),
      99             :     _storeFits(True),
     100             :     _polyOrder(-1), _fitAxis(axis), _nGaussSinglets(0),
     101             :     _nGaussMultiplets(0), _nLorentzSinglets(0), _nPLPCoeffs(0),
     102             :     _nLTPCoeffs(0), _minGoodPoints(1), _nProfiles(0), _nAttempted(0), _nSucceeded(0),
     103             :     _nConverged(0), _nValid(0), _results(Record()), _nonPolyEstimates(),
     104             :     _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 &regionPtr,    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             :     _residual(), _model(), _xUnit(), _centerName(),
     129             :     _centerErrName(), _fwhmName(), _fwhmErrName(),
     130             :     _ampName(), _ampErrName(), _integralName(),
     131             :     _integralErrName(), _plpName(), _plpErrName(), _sigmaName(),
     132             :     _abscissaDivisorForDisplay("1"), _multiFit(False),
     133             :     /* _deleteImageOnDestruct(False),*/ _logResults(True),
     134             :     _isSpectralIndex(False), _createResid(False), _overwrite(overwrite),
     135             :     _storeFits(True),
     136             :     _polyOrder(-1), _fitAxis(axis), _nGaussSinglets(0),
     137             :     _nGaussMultiplets(0), _nLorentzSinglets(0), _nPLPCoeffs(0),
     138             :     _nLTPCoeffs(0), _minGoodPoints(1), _nProfiles(0), _nAttempted(0), _nSucceeded(0),
     139             :     _nConverged(0), _nValid(0), _results(Record()), _nonPolyEstimates(),
     140             :     _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             :         }
     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             :             }
     294             :         }
     295           0 :         _fitallprofiles();
     296           0 :         *_getLog() << logOrigin;
     297             :     }
     298           0 :     catch (const AipsError& x) {
     299           0 :         ThrowCc("Exception during fit: " + x.getMesg());
     300             :     }
     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             : }
     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             :     }
     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             :     }
     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             :         }
     455           0 :         if (mymax != 1) {
     456           0 :             _sigma->put(_sigma->get()/mymax);
     457             :         }
     458             :     }
     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             :         }
     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             :     }
     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             :     }
     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             :     }
     556           0 :     return outputs;
     557             : }
     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             :     }
     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             :     }
     638             :     else {
     639           0 :         convertedVal = world(_fitAxis);
     640             :     }
     641           0 :     return convertedVal;
     642             : }
     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             :     }
     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             :     }
     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             :             std::transform(
     736             :                 _goodPlanes.begin(), _goodPlanes.end(), goodPlanes.begin(),
     737           0 :                 bind2nd(minus<Int>(), imageOff)
     738           0 :             );
     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             :     }
     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             :     }
     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           0 :                 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             :             }
     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             :     }
     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             :             }
    1021           0 :             newEstimates = _fitters(nearest)->getList();
    1022             :         }
    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             :             }
    1091           0 :             default:
    1092           0 :                 continue;
    1093             :             }
    1094             :         }
    1095             :     }
    1096             : }
    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             : }

Generated by: LCOV version 1.16