casa  $Rev:20696$
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
ImageFit1D.h
Go to the documentation of this file.
00001 
00002 //# ImageFit1D.h: Class to fit profiles to vectors from images
00003 //# Copyright (C) 2004
00004 //# Associated Universities, Inc. Washington DC, USA.
00005 //#
00006 //# This library is free software; you can redistribute it and/or modify it
00007 //# under the terms of the GNU Library General Public License as published by
00008 //# the Free Software Foundation; either version 2 of the License, or (at your
00009 //# option) any later version.
00010 //#
00011 //# This library is distributed in the hope that it will be useful, but WITHOUT
00012 //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00013 //# FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Library General Public
00014 //# License for more details.
00015 //#
00016 //# You should have received a copy of the GNU Library General Public License
00017 //# along with this library; if not, write to the Free Software Foundation,
00018 //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
00019 //#
00020 //# Correspondence concerning AIPS++ should be addressed as follows:
00021 //#        Internet email: aips2-request@nrao.edu.
00022 //#        Postal address: AIPS++ Project Office
00023 //#                        National Radio Astronomy Observatory
00024 //#                        520 Edgemont Road
00025 //#                        Charlottesville, VA 22903-2475 USA
00026 //#
00027 //#   $Id: ImageFit1D.h 20229 2008-01-29 15:19:06Z gervandiepen $
00028 
00029 #ifndef IMAGES_IMAGEFIT1D_H
00030 #define IMAGES_IMAGEFIT1D_H
00031 
00032 //# Includes
00033 #include <casa/aips.h>
00034 #include <casa/Arrays/Vector.h>
00035 #include <scimath/Mathematics/NumericTraits.h>
00036 
00037 #include <measures/Measures/MDirection.h>
00038 #include <measures/Measures/MFrequency.h>
00039 #include <measures/Measures/MDoppler.h>
00040 #include <coordinates/Coordinates/CoordinateSystem.h>
00041 #include <components/SpectralComponents/ProfileFit1D.h>
00042 
00043 #include <memory>
00044 
00045 namespace casa {
00046 
00047 class SpectralElement;
00048 class SpectralList;
00049 class ImageRegion;
00050 template<class T> class ImageInterface;
00051 
00052 
00053 // <summary>
00054 // Fit spectral components to a Vector of data from an image
00055 // </summary>
00056 
00057 // <use visibility=export>
00058 
00059 // <reviewed reviewer="" date="" tests="tImageFit1D.cc">
00060 // </reviewed>
00061 
00062 // <prerequisite>
00063 //   <li> <linkto class="SpectralElement">SpectralElement</linkto> 
00064 //   <li> <linkto class="SpectralList">SpectralList</linkto> 
00065 //   <li> <linkto class="SpectralFit">SpectralFit</linkto> 
00066 //   <li> <linkto class="ProfileFit1D">ProfileFit1D</linkto> 
00067 // </prerequisite>
00068 
00069 // <synopsis> 
00070 // Fit lists (held in class SpectralList) of SpectralElements to a Vector of data 
00071 // from the image.  Each SpectralElement can  be one from a variety of types.
00072 // The values of the parameters for each SpectralElement provide the initial 
00073 // starting guesses for the fitting process.  
00074 //
00075 // You specify the domain in which the fit is to be done via the 
00076 // enum AbcissaType.  The CoordinateSystem in the image is used
00077 // to convert the pixel coordinates to the desired abcissa.  
00078 // You can change the units of the CoordinateSystem if you want
00079 // to fit in different units.  If you set an estimate yourself
00080 // (function setElements or addElement) it is the callers responsibility
00081 // that the elements are in the correct abcissa domain.  Function
00082 // setGaussianElements will automatically make an estimate in the 
00083 // correct domain.
00084 //
00085 // Also, a SpectralElement object holds a mask indicating whether 
00086 // a parameter should be held fixed or solved for.   After the 
00087 // fitting is done, a new SpectralList holding SpectralElements with 
00088 // the fitted parameters is created.  
00089 //
00090 // For all the functions that return a status Bool, True is good. If
00091 // False is returned, an error message can be recovered with function
00092 // <src>errorMessage</src>,  You should not proceed if False is returned.
00093 // 
00094 // Exceptions will be thrown if you do not set the Image and axis 
00095 // via the constructor or <src>setImage</src> function.
00096 // </synopsis> 
00097 
00098 // <example>
00099 // <srcblock>
00100 // PagedImage<Float> im("myimage");
00101 // Int axis = 2;
00102 // ImageFit1D<Float> fitter(image, axis);
00103 // IPosition pos(in.ndim(),0);
00104 // fitter.setData(pos, ImageFit1D<Float>::IM_NATIVE);     // Fit in native coordinate space
00105 // fitter.setGaussianElements(3);                      // FIt 3 Gaussians
00106 // if (fitter.fit()) {
00107 //    cerr << fitter.getList() << endl;                // Print result
00108 // }
00109 // 
00110 // </srcblock>
00111 // </example>
00112 
00113 // <todo asof="2004/07/10">
00114 //   <li> Add constraints
00115 // </todo>
00116 
00117 template <class T> class ImageFit1D 
00118 {
00119 public:
00120 
00121     enum AbcissaType {
00122        PIXEL = 0,
00123        IM_NATIVE = 1,
00124        VELOCITY = 2,
00125        N_TYPES};
00126 
00127     // Default Constructor. No image or axis set, so
00128     // <src>fit()</src> cannot be called until other parameters are set,
00129     // via eg <src>setImage</src>.
00130     ImageFit1D();
00131 
00132     // Constructor.  Fitting weights are assumed all unity.
00133     ImageFit1D(const ImageInterface<T>& image, uInt axis=0);
00134 
00135     // Constructor with fitting weights image.  The data and weights images must
00136     // be the same shape.
00137     ImageFit1D(const ImageInterface<T>& image, 
00138                const ImageInterface<T>& weights,
00139                uInt axis=0);
00140 
00141     // Destructor
00142     ~ImageFit1D();
00143 
00144     // Copy constructor.  Uses reference semantics.
00145     ImageFit1D(const ImageFit1D& other);
00146 
00147     // Assignment operator. Uses reference semantics.
00148     ImageFit1D& operator=(const ImageFit1D& other);
00149 
00150     // Set Image(s) and axis
00151     // <group>
00152     void setImage (const ImageInterface<T>& im, const ImageInterface<T>& weights, uInt axis);
00153     void setImage (const ImageInterface<T>& im, uInt axis);
00154     // </group>
00155 
00156     // Set the data to be fit.  All non-profile axes data are averaged.
00157     // For the profile axis, the full spectrum is taken.  The abcissa
00158     // world values are computed when you call these functions.  The domain of the
00159     // abcissa values is controlled by <src>AbcissaType</src> and
00160     // <src>doAbs</src> (absolute coordinates).  The CoordinateSystem in
00161     // the image is used to convert from pixels to world values.
00162     // <group>
00163     Bool setData (
00164         const IPosition& pos, const ImageFit1D<T>::AbcissaType type,
00165         const Bool doAbs=True
00166     );
00167     Bool setData (
00168         const ImageRegion& region, const ImageFit1D<T>::AbcissaType type,
00169         Bool doAbs=True
00170     );
00171     // </group>
00172 
00173     // Set a SpectralList of SpectralElements to fit for.    These elements
00174     // must be in the correct abcissa domain set in function <src>setData</src>.
00175     // You must have already called <src>setData</src> to call this function.
00176     // The SpectralElements in the list hold the
00177     // initial estimates.  They also contain the information about whether
00178     // specific parameters are to be held fixed or allowed to vary in
00179     // the fitting process.
00180     // You can recover the list of elements with function getList.
00181     void setElements (const SpectralList& list) {itsFitter.setElements(list);};
00182 
00183     // Add new SpectralElement(s) to the SpectralList (can be empty)
00184     // of SpectralElements to be fit for.  
00185     // You must have already called <src>setData</src> to call this function.
00186     //<group>
00187     void addElement (const SpectralElement& el) {itsFitter.addElement(el);};
00188     void addElements (const SpectralList& list) {itsFitter.addElements(list);};
00189     // </group>
00190 
00191     // Set a SpectralList of Gaussian SpectralElements to fit for.  
00192     // The initial estimates for the Gaussians will be automatically determined
00193     // in the correct abcissa domain.
00194     // All of the parameters created by this function will be solved for
00195     // by default. You can recover the list of elements with function getList.
00196     // Status is returned, if False, error message can be recovered with <src>errorMessage</src>
00197     Bool setGaussianElements (uInt nGauss);
00198 
00199     // Clear the SpectralList of elements to be fit for
00200     void clearList () {itsFitter.clearList();};
00201 
00202     // Do the fit and return convergence status.  Errors in the fitting
00203     // process will generate an AipsError exception and you should catch
00204     // these yourself.
00205     Bool fit ();
00206 
00207     // Get Chi Squared of fit
00208     Double getChiSquared () const {return itsFitter.getChiSquared();}
00209 
00210     // Get number of iterations for last fit
00211     Double getNumberIterations () const {return itsFitter.getNumberIterations();}
00212 
00213     // Recover the list of elements.  You can get the elements
00214     // as initially estimated (fit=False), or after fitting 
00215     // (fit=True).  In the latter case, the SpectralElements
00216     // hold the parameters and errors of the fit.
00217     const SpectralList& getList (Bool fit=True) const {return itsFitter.getList(fit);};
00218 
00219     // Recover vectors for the estimate, fit and residual.
00220     // If you don't specify which element, all elements are included
00221     // If the Vectors are returned with zero length, it means an error
00222     // condition exists (e.g. asking for fit before you do one). In this
00223     // case an error message can be recovered with function <src>errorMessage</src>.
00224     //<group>
00225     Vector<T> getEstimate (Int which=-1) const;
00226     Vector<T> getFit (Int which=-1) const;
00227     Vector<T> getResidual (Int which=-1, Bool fit=True)  const;
00228     //</group>
00229 
00230     // Get Total Mask (data and range mask)
00231     Vector<Bool> getTotalMask () const {return itsFitter.getTotalMask();};
00232 
00233     // did the fit succeed? should only be called after fit().
00234     Bool succeeded() const;
00235 
00236     // did the fit converge? should only be called after fit().
00237     Bool converged() const;
00238 
00239     // Recover the error message
00240     String errorMessage () const {return itsError;};
00241 
00242     // Helper function.  Sets up the CoordinateSystem to reflect the choice of
00243     // abcissa unit and the doppler (if the axis is spectral).
00244     static Bool setAbcissaState (String& errMsg, ImageFit1D<T>::AbcissaType& type,
00245                                  CoordinateSystem& cSys, const String& xUnit,
00246                                  const String& doppler, uInt pixelAxis);
00247 
00248 
00249     // flag the solution as invalid based on external criteria.
00250     void invalidate();
00251 
00252     // is the solution valid? If False, some external logic has
00253     // called invalidate()
00254     Bool isValid() const;
00255 
00256 private:
00257    std::auto_ptr<ImageInterface<T> > itsImagePtr;
00258    std::auto_ptr<ImageInterface<T> > itsWeightPtr;
00259    uInt itsAxis;
00260 
00261 // In the future I will be able to template the fitter on T. For now
00262 // it must be Double.
00263 
00264    typedef typename NumericTraits<T>::PrecisionType FitterType;
00265    ProfileFit1D<FitterType> itsFitter;
00266    CoordinateSystem itsCS;
00267    mutable String itsError;                // Error message
00268    Bool _converged, _success, _isValid;
00269 // Functions
00270    
00271    void check() const;
00272    void checkType() const;
00273    void copy (const ImageFit1D<T>& other);
00274    Bool makeAbcissa (Vector<Double>& x, ImageFit1D<T>::AbcissaType type, Bool doAbs);
00275    void setWeightsImage (const ImageInterface<T>& im);
00276 
00277    // reset the fitter, for example if we've done a fit and want to move
00278    // to the next position in the image
00279    void _resetFitter();
00280 
00281 };
00282 
00283 } //#End casa namespace
00284 
00285 #ifndef CASACORE_NO_AUTO_TEMPLATES
00286 #include <imageanalysis/ImageAnalysis/ImageFit1D.tcc>
00287 #endif //# CASACORE_NO_AUTO_TEMPLATES
00288 #endif