Line data Source code
1 : //# tSubImage.cc: Test program for class SubImage
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: tSubImage.cc 20567 2009-04-09 23:12:39Z gervandiepen $
27 :
28 : #ifndef IMAGES_IMAGEPROFILEFITTER_H
29 : #define IMAGES_IMAGEPROFILEFITTER_H
30 :
31 : #include <imageanalysis/ImageAnalysis/ImageTask.h>
32 :
33 : #include <components/SpectralComponents/GaussianMultipletSpectralElement.h>
34 : #include <imageanalysis/ImageAnalysis/ImageFit1D.h>
35 : #include <casacore/images/Images/TempImage.h>
36 :
37 : #include <casacore/casa/namespace.h>
38 :
39 : namespace casa {
40 :
41 : class ProfileFitResults;
42 :
43 : class ImageProfileFitter : public ImageTask<casacore::Float> {
44 : // <summary>
45 : // Top level interface for one-dimensional profile fits.
46 : // </summary>
47 :
48 : // <reviewed reviewer="" date="" tests="" demos="">
49 : // </reviewed>
50 :
51 : // <prerequisite>
52 : // </prerequisite>
53 :
54 : // <etymology>
55 : // Fits one-dimensional profiles.
56 : // </etymology>
57 :
58 : // <synopsis>
59 : // Top level interface for one-dimensional profile fits.
60 : // </synopsis>
61 :
62 : // <example>
63 : // <srcblock>
64 : // ImageProfileFitter fitter(...)
65 : // fitter.fit()
66 : // </srcblock>
67 : // </example>
68 :
69 : public:
70 : // constructor appropriate for API calls.
71 : // Parameters:
72 : // <src>image</src> - the input image in which to fit the models
73 : // <src>box</src> - A 2-D rectangular box in direction space in which to use pixels for the fitting, eg box=100,120,200,230
74 : // In cases where both box and region are specified, box, not region, is used.
75 : // <src>region</src> - Named region to use for fitting. "" => Don't use a named region
76 : // <src>regPtr</src> - Pointer to a region record. 0 => don't use a region record.
77 : // <src>chans</src> - Zero-based channel range on which to do the fit.
78 : // <src>stokes</src> - casacore::Stokes plane on which to do the fit. Only a single casacore::Stokes parameter can be
79 : // specified.
80 : // Only a maximum of one of region, regionPtr, or box/stokes/chans should be specified.
81 : // <src>mask</src> - Mask (as LEL) to use as a way to specify which pixels to use </src>
82 : // <src>axis</src> - axis along which to do the fits. If <0, use spectral axis, and if no spectral
83 : // axis, use zeroth axis.
84 : // <src>ngauss</src> number of single gaussians to fit.
85 : // <src>estimatesFilename</src> file containing initial estimates for single gaussians.
86 : // <src>spectralList</src> spectral list containing initial estimates of single gaussians. Do
87 : // not put a polynomial in here; set that with setPolyOrder().
88 :
89 : // Fit only Gaussian singlets and an optional polynomial. In this case, the
90 : // code guestimates initial estimates for the specified number of Gaussian
91 : // singlets. If you only wish to fit a polynomial, you must use this
92 : // constructor and you must set <src>ngauss</src> to zero. After construction,
93 : // you must call setPolyOrder().
94 : ImageProfileFitter(
95 : const SPCIIF image, const casacore::String& region,
96 : const casacore::Record *const ®ionPtr, const casacore::String& box,
97 : const casacore::String& chans, const casacore::String& stokes, const casacore::String& mask,
98 : const casacore::Int axis, const casacore::uInt ngauss, casacore::Bool overwrite=false
99 : );
100 :
101 : // Fit only Gaussian singlets and an optional polynomial. In this case, the number
102 : // of Gaussian singlets is deduced from the specified estimates file.
103 : ImageProfileFitter(
104 : const SPCIIF image, const casacore::String& region,
105 : const casacore::Record *const ®ionPtr, const casacore::String& box,
106 : const casacore::String& chans, const casacore::String& stokes, const casacore::String& mask,
107 : const casacore::Int axis, const casacore::String& estimatesFilename,
108 : casacore::Bool overwrite=false
109 : );
110 :
111 : // Fit any permitted combination of spectral components and an optional polynomial.
112 : // All components to be fit (except a possible polynomial) must be represented
113 : // with initial estimates in <src>spectralList</src>.
114 : ImageProfileFitter(
115 : const SPCIIF image, const casacore::String& region,
116 : const casacore::Record *const ®ionPtr, const casacore::String& box,
117 : const casacore::String& chans, const casacore::String& stokes, const casacore::String& mask,
118 : const casacore::Int axis, const SpectralList& spectralList, casacore::Bool overwrite=false
119 : );
120 :
121 : // destructor
122 : ~ImageProfileFitter();
123 :
124 : // Do the fit. If doDetailedResults is false, an empty casacore::Record is returned.
125 : casacore::Record fit(casacore::Bool doDetailedResults=true);
126 :
127 : // get the fit results. If fit was run with doDetailedResults=false, an empty casacore::Record is returned
128 : casacore::Record getResults() const;
129 :
130 32 : inline casacore::String getClass() const { return _class; };
131 :
132 : // set the order of a polynomial to be simultaneously fit.
133 : void setPolyOrder(casacore::Int p);
134 :
135 : // set whether to do a pixel by pixel fit.
136 : inline void setDoMultiFit(const casacore::Bool m) { _multiFit = m; }
137 :
138 : // set if results should be written to the logger
139 : inline void setLogResults(const casacore::Bool logResults) { _logResults = logResults; }
140 :
141 : // set minimum number of good points required to attempt a fit
142 : inline void setMinGoodPoints(const casacore::uInt mgp) {
143 : ThrowIf(mgp == 0, "Number of good points has to be > 0");
144 : _minGoodPoints = mgp;
145 : }
146 :
147 : // <group>
148 : // Solution images. Only written if _multifit is true
149 : // model image name
150 : inline void setModel(const casacore::String& model) { _model = model; }
151 : // residual image name
152 : inline void setResidual(const casacore::String& residual) { _residual = residual; }
153 : // gaussian amplitude image name
154 : inline void setAmpName(const casacore::String& s) { _ampName = s; }
155 : // gaussian amplitude error image name
156 : inline void setAmpErrName(const casacore::String& s) { _ampErrName = s; }
157 : // gaussian center image name
158 : inline void setCenterName(const casacore::String& s) { _centerName = s; }
159 : // gaussian center error image name
160 : inline void setCenterErrName(const casacore::String& s) { _centerErrName = s; }
161 : // gaussian fwhm image name
162 : inline void setFWHMName(const casacore::String& s) { _fwhmName = s; }
163 : // gaussian fwhm error image name
164 : inline void setFWHMErrName(const casacore::String& s) { _fwhmErrName = s; }
165 : // gaussian integral image name
166 : inline void setIntegralName(const casacore::String& s) { _integralName = s; }
167 : // gaussian integral error image name
168 : inline void setIntegralErrName(const casacore::String& s) { _integralErrName = s; }
169 : // </group>
170 :
171 : // set the name of the power logarithmic polynomial image.
172 : inline void setPLPName(const casacore::String& s) { _plpName = s; }
173 :
174 : // set the name of the power logarithmic polynomial image.
175 : inline void setPLPErrName(const casacore::String& s) { _plpErrName = s; }
176 :
177 : // set the name of the logarithmic transformed polynomial image.
178 : inline void setLTPName(const casacore::String& s) { _ltpName = s; }
179 :
180 : // set the name of the logarithmic transformed polynomial image.
181 : inline void setLTPErrName(const casacore::String& s) { _ltpErrName = s; }
182 :
183 : // set the range over which PFC amplitude solutions are valid
184 : void setGoodAmpRange(const casacore::Double min, const casacore::Double max);
185 :
186 : // set the range over which PFC center solutions are valid
187 : void setGoodCenterRange(const casacore::Double min, const casacore::Double max);
188 :
189 : // set the range over which PFC FWHM solutions are valid
190 : void setGoodFWHMRange(const casacore::Double min, const casacore::Double max);
191 :
192 : // <group>
193 : // set standard deviation image
194 : void setSigma(const casacore::Array<casacore::Float>& sigma);
195 :
196 : void setSigma(const casacore::ImageInterface<casacore::Float> *const &sigma);
197 :
198 : inline void setOutputSigmaImage(const casacore::String& s) { _sigmaName = s; }
199 : // </group>
200 :
201 : const casacore::Array<std::shared_ptr<ProfileFitResults> >& getFitters() const;
202 :
203 : //Converts a pixel value into a world value either in velocity, wavelength, or
204 : //frequency units. If the tabular index >= 0, it uses the tabular index for conversion
205 : //with the specified casacore::MFrequency type, otherwise, it uses the spectral axis for
206 : //conversion.
207 : casacore::Double getWorldValue(
208 : casacore::Double pixelVal, const casacore::IPosition& imPos, const casacore::String& units,
209 : casacore::Bool velocity, casacore::Bool wavelength, casacore::Int tabularIndex=-1,
210 : casacore::MFrequency::Types type=casacore::MFrequency::DEFAULT
211 : ) const;
212 :
213 : void setAbscissaDivisor(casacore::Double d);
214 :
215 : void setAbscissaDivisor(const casacore::Quantity& q);
216 :
217 : // for backward compatibility with ia.continuumsub. If no residual
218 : // image has been provided, a casacore::TempImage is created.
219 : void createResidualImage(casacore::Bool b) { _createResid = b; }
220 :
221 : SPIIF getResidual() const {
222 : return _residImage;
223 : }
224 :
225 : // set the planes along the fit axis that are considered good for performing
226 : // the fits. An empty vector implies that all planes are considered "good".
227 : void setGoodPlanes(const std::set<casacore::uInt>& planes) { _goodPlanes = planes; }
228 :
229 : protected:
230 :
231 96 : inline CasacRegionManager::StokesControl _getStokesControl() const {
232 96 : return CasacRegionManager::USE_FIRST_STOKES;
233 : }
234 :
235 96 : inline std::vector<casacore::Coordinate::Type> _getNecessaryCoordinates() const {
236 96 : return std::vector<casacore::Coordinate::Type>(0);
237 : }
238 :
239 96 : casacore::Bool _hasLogfileSupport() const { return true; }
240 :
241 96 : inline casacore::Bool _supportsMultipleRegions() const {return true;}
242 :
243 : private:
244 : casacore::String _residual, _model, _xUnit,
245 : _centerName, _centerErrName, _fwhmName,
246 : _fwhmErrName, _ampName, _ampErrName,
247 : _integralName, _integralErrName, _plpName, _plpErrName,
248 : _ltpName, _ltpErrName, _sigmaName, _abscissaDivisorForDisplay;
249 : casacore::Bool _multiFit, _logResults, _isSpectralIndex,
250 : _createResid, _overwrite, _storeFits;
251 : casacore::Int _polyOrder, _fitAxis;
252 : casacore::uInt _nGaussSinglets, _nGaussMultiplets, _nLorentzSinglets,
253 : _nPLPCoeffs, _nLTPCoeffs, _minGoodPoints, _nProfiles, _nAttempted,
254 : _nSucceeded, _nConverged, _nValid;
255 : casacore::Array<std::shared_ptr<ProfileFitResults> > _fitters;
256 : // subimage contains the region of the original image
257 : // on which the fit is performed.
258 : std::shared_ptr<const casacore::SubImage<casacore::Float> > _subImage;
259 : casacore::Record _results;
260 : SpectralList _nonPolyEstimates;
261 : std::unique_ptr<std::pair<casacore::Double, casacore::Double> > _goodAmpRange, _goodCenterRange, _goodFWHMRange;
262 : casacore::Matrix<casacore::String> _worldCoords;
263 : std::shared_ptr<casacore::TempImage<casacore::Float> > _sigma;
264 : casacore::Double _abscissaDivisor;
265 : SPIIF _modelImage, _residImage;
266 : // planes along _fitAxis to use in fits, empty => use all planes
267 : // originally used to support continuum subtraction
268 : std::set<casacore::uInt> _goodPlanes;
269 :
270 : const static casacore::String _class;
271 :
272 : mutable casacore::Bool _haveWarnedAboutGuessingGaussians = false;
273 :
274 : std::vector<OutputDestinationChecker::OutputStruct> _getOutputStruct();
275 :
276 : void _checkNGaussAndPolyOrder() const;
277 :
278 : void _finishConstruction();
279 :
280 : // moved from ImageAnalysis
281 : void _fitallprofiles();
282 :
283 : // Fit all profiles in image. The output images must be already
284 : // created; if the pointer is 0, that image won't be filled.
285 : // The mask from the input image is transferred to the output
286 : // images. If the weights image is pointer is non-zero, the
287 : // values from it will be used to weight the data points in the
288 : // fit. You can fit some combination of gaussians and a polynomial
289 : // (-1 means no polynomial). Initial estimates are not required.
290 : // Fits are done in image space to provide astronomer friendly results,
291 : // but pixel space is better for the fitter when fitting polynomials.
292 : // Thus, atm, callers should be aware that fitting polynomials may
293 : // fail even when the data lie exactly on a polynomial curve.
294 : // This will probably be fixed in the future by doing the fits
295 : // in pixel space here and requiring the caller to deal with converting
296 : // to something astronomer friendly if it so desires.
297 :
298 : void _fitProfiles(const casacore::Bool showProgress);
299 :
300 : //casacore::Bool _inVelocitySpace() const;
301 :
302 : void _flagFitterIfNecessary(ImageFit1D<casacore::Float>& fitter) const;
303 :
304 : casacore::Bool _isPCFSolutionOK(const PCFSpectralElement *const &pcf) const;
305 :
306 : void _loopOverFits(
307 : SPCIIF fitData, casacore::Bool showProgress,
308 : std::shared_ptr<casacore::ProgressMeter> progressMeter, casacore::Bool checkMinPts,
309 : const casacore::Array<casacore::Bool>& fitMask, ImageFit1D<casacore::Float>::AbcissaType abcissaType,
310 : const casacore::IPosition& fitterShape, const casacore::IPosition& sliceShape,
311 : const std::set<casacore::uInt> goodPlanes
312 : );
313 :
314 : void _setAbscissaDivisorIfNecessary(const casacore::Vector<casacore::Double>& abscissaValues);
315 :
316 : casacore::Bool _setFitterElements(
317 : ImageFit1D<casacore::Float>& fitter, SpectralList& newEstimates,
318 : const std::unique_ptr<const PolynomialSpectralElement>& polyEl,
319 : const std::vector<casacore::IPosition>& goodPos,
320 : const casacore::IPosition& fitterShape, const casacore::IPosition& curPos,
321 : casacore::uInt nOrigComps
322 : ) const;
323 :
324 : void _updateModelAndResidual(
325 : casacore::Bool fitOK, const ImageFit1D<casacore::Float>& fitter,
326 : const casacore::IPosition& sliceShape, const casacore::IPosition& curPos,
327 : casacore::Lattice<casacore::Bool>* const &pFitMask, casacore::Lattice<casacore::Bool>* const &pResidMask
328 : ) const;
329 :
330 : // only implemented for the simplest cases so far
331 : casacore::uInt _nUnknowns() const;
332 :
333 : };
334 : }
335 :
336 : #endif
|
