Line data Source code
1 :
2 : //# ImageFit1D.h: Class to fit profiles to vectors from images
3 : //# Copyright (C) 2004
4 : //# Associated Universities, Inc. Washington DC, USA.
5 : //#
6 : //# This library is free software; you can redistribute it and/or modify it
7 : //# under the terms of the GNU Library General Public License as published by
8 : //# the Free Software Foundation; either version 2 of the License, or (at your
9 : //# option) any later version.
10 : //#
11 : //# This library is distributed in the hope that it will be useful, but WITHOUT
12 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
13 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
14 : //# License for more details.
15 : //#
16 : //# You should have received a copy of the GNU Library General Public License
17 : //# along with this library; if not, write to the Free Software Foundation,
18 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
19 : //#
20 : //# Correspondence concerning AIPS++ should be addressed as follows:
21 : //# Internet email: aips2-request@nrao.edu.
22 : //# Postal address: AIPS++ Project Office
23 : //# National Radio Astronomy Observatory
24 : //# 520 Edgemont Road
25 : //# Charlottesville, VA 22903-2475 USA
26 : //#
27 : //# $Id: ImageFit1D.h 20229 2008-01-29 15:19:06Z gervandiepen $
28 :
29 : #ifndef IMAGES_IMAGEFIT1D_H
30 : #define IMAGES_IMAGEFIT1D_H
31 :
32 : //# Includes
33 : #include <casacore/casa/aips.h>
34 : #include <casacore/casa/Arrays/Vector.h>
35 : #include <casacore/scimath/Mathematics/NumericTraits.h>
36 : #include <casacore/measures/Measures/MDirection.h>
37 : #include <casacore/measures/Measures/MFrequency.h>
38 : #include <casacore/measures/Measures/MDoppler.h>
39 : #include <casacore/coordinates/Coordinates/CoordinateSystem.h>
40 : #include <components/SpectralComponents/ProfileFit1D.h>
41 :
42 : #include <memory>
43 :
44 : namespace casacore{
45 :
46 : class ImageRegion;
47 : template<class T> class ImageInterface;
48 : }
49 :
50 : namespace casa {
51 :
52 : class SpectralElement;
53 : class SpectralList;
54 :
55 :
56 : // <summary>
57 : // Fit spectral components to a casacore::Vector of data from an image
58 : // </summary>
59 :
60 : // <use visibility=export>
61 :
62 : // <reviewed reviewer="" date="" tests="tImageFit1D.cc">
63 : // </reviewed>
64 :
65 : // <prerequisite>
66 : // <li> <linkto class="SpectralElement">SpectralElement</linkto>
67 : // <li> <linkto class="SpectralList">SpectralList</linkto>
68 : // <li> <linkto class="SpectralFit">SpectralFit</linkto>
69 : // <li> <linkto class="ProfileFit1D">ProfileFit1D</linkto>
70 : // </prerequisite>
71 :
72 : // <synopsis>
73 : // Fit lists (held in class SpectralList) of SpectralElements to a casacore::Vector of data
74 : // from the image. Each SpectralElement can be one from a variety of types.
75 : // The values of the parameters for each SpectralElement provide the initial
76 : // starting guesses for the fitting process.
77 : //
78 : // You specify the domain in which the fit is to be done via the
79 : // enum AbcissaType. The casacore::CoordinateSystem in the image is used
80 : // to convert the pixel coordinates to the desired abcissa.
81 : // You can change the units of the casacore::CoordinateSystem if you want
82 : // to fit in different units. If you set an estimate yourself
83 : // (function setElements or addElement) it is the callers responsibility
84 : // that the elements are in the correct abcissa domain. Function
85 : // setGaussianElements will automatically make an estimate in the
86 : // correct domain.
87 : //
88 : // Also, a SpectralElement object holds a mask indicating whether
89 : // a parameter should be held fixed or solved for. After the
90 : // fitting is done, a new SpectralList holding SpectralElements with
91 : // the fitted parameters is created.
92 : //
93 : // For all the functions that return a status casacore::Bool, true is good. If
94 : // false is returned, an error message can be recovered with function
95 : // <src>errorMessage</src>, You should not proceed if false is returned.
96 : //
97 : // Exceptions will be thrown if you do not set the Image and axis
98 : // via the constructor or <src>setImage</src> function.
99 : // </synopsis>
100 :
101 : // <example>
102 : // <srcblock>
103 : // casacore::PagedImage<casacore::Float> im("myimage");
104 : // casacore::Int axis = 2;
105 : // ImageFit1D<casacore::Float> fitter(image, axis);
106 : // casacore::IPosition pos(in.ndim(),0);
107 : // fitter.setData(pos, ImageFit1D<casacore::Float>::IM_NATIVE); // Fit in native coordinate space
108 : // fitter.setGaussianElements(3); // FIt 3 Gaussians
109 : // if (fitter.fit()) {
110 : // cerr << fitter.getList() << endl; // Print result
111 : // }
112 : //
113 : // </srcblock>
114 : // </example>
115 :
116 : // <todo asof="2004/07/10">
117 : // <li> Add constraints
118 : // </todo>
119 :
120 :
121 : template <class T> class ImageFit1D {
122 : public:
123 :
124 : using FitterType = typename casacore::NumericTraits<T>::PrecisionType;
125 :
126 : enum AbcissaType {
127 : PIXEL = 0,
128 : IM_NATIVE = 1,
129 : VELOCITY = 2,
130 : N_TYPES};
131 :
132 :
133 : // Constructor. Fitting weights are assumed all unity.
134 : ImageFit1D(std::shared_ptr<const casacore::ImageInterface<T> > image, casacore::uInt axis=0);
135 :
136 : // Constructor with fitting weights image. The data and weights images must
137 : // be the same shape.
138 : ImageFit1D(
139 : std::shared_ptr<const casacore::ImageInterface<T> > image,
140 : std::shared_ptr<const casacore::ImageInterface<T> > weights, casacore::uInt axis=0
141 : );
142 :
143 : // Destructor
144 : ~ImageFit1D();
145 :
146 : // Copy constructor. Uses reference semantics.
147 : ImageFit1D(const ImageFit1D& other);
148 :
149 : // Assignment operator. Uses reference semantics.
150 : ImageFit1D& operator=(const ImageFit1D& other);
151 :
152 :
153 : // Set the data to be fit. All non-profile axes data are averaged.
154 : // For the profile axis, the full spectrum is taken. The abscissa
155 : // world values are computed when you call these functions unless they
156 : // have been set previously by a call to setAbscissa() in which case
157 : // the values that were passed to that method are used. Use the first
158 : // form of setData() in this case. The domain of the
159 : // abscissa values is controlled by <src>AbcissaType</src> and
160 : // <src>doAbs</src> (absolute coordinates). The casacore::CoordinateSystem in
161 : // the image is used to convert from pixels to world values.
162 : // If <src>type</src>=IN_NATIVE and <src>abscissaDivisor</src> is not null,
163 : // the world abscissa values will be divided by the value pointed to by
164 : // <src>abscissaDivisor</src>. This mitigates having very large or very small
165 : // abscissa values when fitting. If xfunc and/or yfunc is not NULL, the x and/or
166 : // y values are fed to the specified function and the resultant values are what
167 : // are used for the x and/or y values in the fit. If xfunc is not NULL and
168 : // setAbscissa values has been called prior, no abscissa value transformation occurs.
169 : // Thus if you want to apply a function to the abscissa values, the caller should
170 : // pass the result of that function into setAbscissaValues.
171 : // <group>
172 : void setData (
173 : const casacore::IPosition& pos, /*const ImageFit1D<T>::AbcissaType type,
174 : const casacore::Bool doAbs=true, const casacore::Double* const &abscissaDivisor=0,
175 : casacore::Array<casacore::Double> (*xfunc)(const casacore::Array<casacore::Double>&)=0, */
176 : casacore::Array<FitterType> (*yfunc)(const casacore::Array<FitterType>&)=0
177 : );
178 :
179 : void setData (
180 : const casacore::IPosition& pos, const ImageFit1D<T>::AbcissaType type,
181 : const casacore::Bool doAbs=true, const casacore::Double* const &abscissaDivisor=0,
182 : casacore::Array<casacore::Double> (*xfunc)(const casacore::Array<casacore::Double>&)=0,
183 : casacore::Array<FitterType> (*yfunc)(const casacore::Array<FitterType>&)=0
184 : );
185 :
186 : /*
187 : casacore::Bool setData (
188 : const casacore::ImageRegion& region, const ImageFit1D<T>::AbcissaType type,
189 : casacore::Bool doAbs=true
190 : );
191 : */
192 : // </group>
193 :
194 : // Set a SpectralList of SpectralElements to fit for. These elements
195 : // must be in the correct abcissa domain set in function <src>setData</src>.
196 : // You must have already called <src>setData</src> to call this function.
197 : // The SpectralElements in the list hold the
198 : // initial estimates. They also contain the information about whether
199 : // specific parameters are to be held fixed or allowed to vary in
200 : // the fitting process.
201 : // You can recover the list of elements with function getList.
202 544 : void setElements (const SpectralList& list) {_fitter.setElements(list);};
203 :
204 : // Add new SpectralElement(s) to the SpectralList (can be empty)
205 : // of SpectralElements to be fit for.
206 : // You must have already called <src>setData</src> to call this function.
207 : //<group>
208 529010 : void addElement (const SpectralElement& el) {_fitter.addElement(el);};
209 : void addElements (const SpectralList& list) {_fitter.addElements(list);};
210 : // </group>
211 :
212 : // Set a SpectralList of Gaussian SpectralElements to fit for.
213 : // The initial estimates for the Gaussians will be automatically determined
214 : // in the correct abcissa domain.
215 : // All of the parameters created by this function will be solved for
216 : // by default. You can recover the list of elements with function getList.
217 : // Status is returned, if false, error message can be recovered with <src>errorMessage</src>
218 : void setGaussianElements (casacore::uInt nGauss);
219 :
220 : // Clear the SpectralList of elements to be fit for
221 545684 : void clearList () {_fitter.clearList();};
222 :
223 : // Do the fit and return convergence status. Errors in the fitting
224 : // process will generate an casacore::AipsError exception and you should catch
225 : // these yourself.
226 : casacore::Bool fit ();
227 :
228 : // Get Chi Squared of fit
229 : casacore::Double getChiSquared () const {return _fitter.getChiSquared();}
230 :
231 : // Get number of iterations for last fit
232 16755 : casacore::Double getNumberIterations () const {return _fitter.getNumberIterations();}
233 :
234 : // Recover the list of elements. You can get the elements
235 : // as initially estimated (fit=false), or after fitting
236 : // (fit=true). In the latter case, the SpectralElements
237 : // hold the parameters and errors of the fit.
238 591300 : const SpectralList& getList (casacore::Bool fit=true) const {return _fitter.getList(fit);};
239 :
240 : // Recover vectors for the estimate, fit and residual.
241 : // If you don't specify which element, all elements are included
242 : // If the Vectors are returned with zero length, it means an error
243 : // condition exists (e.g. asking for fit before you do one). In this
244 : // case an error message can be recovered with function <src>errorMessage</src>.
245 : //<group>
246 : casacore::Vector<T> getEstimate (casacore::Int which=-1) const;
247 : casacore::Vector<T> getFit (casacore::Int which=-1) const;
248 : casacore::Vector<T> getResidual (casacore::Int which=-1, casacore::Bool fit=true) const;
249 : //</group>
250 :
251 4641 : casacore::Bool setXMask(const std::set<casacore::uInt>& indices, casacore::Bool specifiedPixelsAreGood) {
252 4641 : return _fitter.setXMask(indices, specifiedPixelsAreGood);
253 : }
254 :
255 : // get data mask
256 529098 : casacore::Vector<casacore::Bool> getDataMask () const {return _fitter.getDataMask();};
257 :
258 : // Get Total Mask (data and range mask)
259 : casacore::Vector<casacore::Bool> getTotalMask () const {return _fitter.getTotalMask();};
260 :
261 : // did the fit succeed? should only be called after fit().
262 : casacore::Bool succeeded() const;
263 :
264 : // did the fit converge? should only be called after fit().
265 : casacore::Bool converged() const;
266 :
267 : // Helper function. Sets up the casacore::CoordinateSystem to reflect the choice of
268 : // abcissa unit and the doppler (if the axis is spectral).
269 : static casacore::Bool setAbcissaState (casacore::String& errMsg, ImageFit1D<T>::AbcissaType& type,
270 : casacore::CoordinateSystem& cSys, const casacore::String& xUnit,
271 : const casacore::String& doppler, casacore::uInt pixelAxis);
272 :
273 :
274 : // flag the solution as invalid based on external criteria.
275 : void invalidate();
276 :
277 : // is the solution valid? If false, some external logic has
278 : // called invalidate()
279 : casacore::Bool isValid() const;
280 :
281 : // Set the abscissa values prior to running setData. If this is done, then
282 : // the abscissa values will not be recomputed when setData is called.
283 : // This can imporove performance if, for example, you are looping over several fitters for
284 : // which you know the abscissa values do not change.
285 86 : void setAbscissa(const casacore::Vector<casacore::Double>& x) { _x.assign(x); }
286 :
287 : // make the abscissa values, <src>x</src>. If <src>type</src>=IN_NATIVE
288 : // and <src>abscissaDivisor is not null, then divide the native values
289 : // by the value pointed to by <src>abscissaDivisor</src> in making the abscissa
290 : // values.
291 : casacore::Vector<casacore::Double> makeAbscissa (
292 : ImageFit1D<T>::AbcissaType type,
293 : casacore::Bool doAbs, const casacore::Double* const &abscissaDivisor
294 : );
295 : private:
296 : std::shared_ptr<const casacore::ImageInterface<T> > _image, _weights;
297 : casacore::uInt _axis;
298 :
299 : // In the future I will be able to template the fitter on T. For now
300 : // it must be Double.
301 :
302 : ProfileFit1D<FitterType> _fitter;
303 : casacore::Bool _converged, _success, _isValid;
304 : casacore::Vector<casacore::Double> _x, _unityWeights, _weightSlice;
305 : casacore::IPosition _sliceShape;
306 :
307 : // Disallow default constructor
308 : ImageFit1D() {}
309 :
310 : void check() const;
311 : void checkType() const;
312 : void _construct();
313 : void copy (const ImageFit1D<T>& other);
314 :
315 : //void setWeightsImage (const casacore::ImageInterface<T>& im);
316 :
317 : // reset the fitter, for example if we've done a fit and want to move
318 : // to the next position in the image
319 : void _resetFitter();
320 :
321 :
322 : // Set Image(s) and axis
323 : // <group>
324 : // void setImage (const casacore::ImageInterface<T>& im, const casacore::ImageInterface<T>& weights, casacore::uInt pixelAxis);
325 : // void setImage (const casacore::ImageInterface<T>& im, casacore::uInt pixelAxis);
326 : // </group>
327 :
328 : };
329 :
330 : } //#End casa namespace
331 :
332 : #ifndef CASACORE_NO_AUTO_TEMPLATES
333 : #include <imageanalysis/ImageAnalysis/ImageFit1D.tcc>
334 : #endif //# CASACORE_NO_AUTO_TEMPLATES
335 : #endif
|
