ImageDecomposer.h

Go to the documentation of this file.
00001 //# ImageDecomposer.h: decompose images into components
00002 //# Copyright (C) 1994,1995,1996,1997,1998,1999,2000,2001,2002
00003 //# Associated Universities, Inc. Washington DC, USA.
00004 //#
00005 //# This library is free software; you can redistribute it and/or modify it
00006 //# under the terms of the GNU Library General Public License as published by
00007 //# the Free Software Foundation; either version 2 of the License, or (at your
00008 //# option) any later version.
00009 //#
00010 //# This library is distributed in the hope that it will be useful, but WITHOUT
00011 //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00012 //# FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Library General Public
00013 //# License for more details.
00014 //#
00015 //# You should have received a copy of the GNU Library General Public License
00016 //# along with this library; if not, write to the Free Software Foundation,
00017 //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
00018 //#
00019 //# Correspondence concerning AIPS++ should be addressed as follows:
00020 //#        Internet email: aips2-request@nrao.edu.
00021 //#        Postal address: AIPS++ Project Office
00022 //#                        National Radio Astronomy Observatory
00023 //#                        520 Edgemont Road
00024 //#                        Charlottesville, VA 22903-2475 USA
00025 //#
00026 //# $Id: ImageDecomposer.h 20739 2009-09-29 01:15:15Z Malte.Marquarding $
00027 
00028 #ifndef IMAGES_IMAGEDECOMPOSER_H
00029 #define IMAGES_IMAGEDECOMPOSER_H
00030 
00031 #include <casa/iostream.h>
00032 #include <casa/math.h>
00033 #include <casa/aips.h>
00034 #include <casa/Arrays/Vector.h>
00035 #include <casa/Arrays/Matrix.h>
00036 #include <casa/Containers/Block.h>
00037 #include <scimath/Functionals/Function1D.h>
00038 #include <lattices/Lattices/TempLattice.h>
00039 #include <lattices/Lattices/SubLattice.h>
00040 #include <images/Images/ImageInterface.h>
00041 
00042 
00043 namespace casa { //# NAMESPACE CASA - BEGIN
00044 
00045 // <summary>
00046 // A tool to separate a complex image into individual components.
00047 // </summary>
00048 //
00049 // <use visibility=export>
00050 //
00051 // <reviewed reviewer="" date="" tests="tImageDecomposer.cc">
00052 // </reviewed>
00053 //
00054 // <prerequisite>
00055 //   <li> <linkto class=CoordinateSystem>CoordinateSystem</linkto>   
00056 //   <li> <linkto class=ImageInterface>ImageInterface</linkto>
00057 // </prerequisite>
00058 
00059 // <etymology>
00060 // It takes an image, and separates it into components.
00061 // </etymology>
00062 
00063 // <synopsis>
00064 // ImageDecomposer is an image decomposition tool that performs several tasks,
00065 // with the end result being that a strongly blended image is separated into
00066 // components - both in the sense that it determines the parameters for each
00067 // component (assuming a Gaussian model) and that it physically assigns each
00068 // pixel in the image to an individual object.  The products of these two
00069 // operations are called the component list and the component map, 
00070 // respectively.  The fitting process (which determines the component list) and
00071 // the pixel-decomposition process (which determines the component map) are
00072 // designed to work cooperatively to increase the efficiency and accuracy of
00073 // both, though each can operate without the other if necessary.
00074 // 
00075 // The algorithm between the decomposition is based on the function clfind
00076 // described in Williams et al 1994, which uses a contouring procedure whereby
00077 // a closed contour designates a separate component.  The program first 
00078 // separates the image into clearly distint 'regions' of blended emission, then
00079 // contours each region to determine the areas constituting each component and
00080 // passes this information on to the fitter, which determines the component 
00081 // list.  
00082 //
00083 // The software is compatible with 2 and 3 dimensional images, but is not
00084 // yet structured for higher dimensions.
00085 // </synopsis>
00086 
00087 // <example>
00088 // <srcblock>
00089 //  TempImage<Double> image;
00090 //  //(populate the image with data: see dImageDecomposer.cc)
00091 //  ImageDecomposer<Double> id(image);
00092 //  id.setDeblendOptions(0.3, 8);
00093 //  id.setFitOptions(0.4);
00094 //  id.decomposeImage();
00095 //  id.display();
00096 //  id.printComponents();
00097 // </srcblock>
00098 // </example>
00099 //   
00100 // <motivation>
00101 // </motivation>
00102 //
00103 // <note>
00104 // </note>
00105 //
00106 // <todo asof="2002/07/23">
00107 //   <li> Generalize dimensionality 
00108 //   <li> Use Lattice iterators in place of IPosition loops wherever possible
00109 //   <li> Speed up fitting by not sending every region pixel to the fitter
00110 //   <li> Send the completed componentmap to the user as an AIPS++ (Int?) Image
00111 //   <li> Return a ComponentList instead of a Matrix
00112 //   <li> Enable custom contouring at user level
00113 //   <li> Add progress meter
00114 //   <li> Numerous other improvements to make are documented in the code
00115 // </todo>
00116   
00117 
00118 
00119 template <class T> class ImageDecomposer {
00120 
00121 public:
00122 
00123 // 'Special' flag values for pixels in the component map.  An indeterminate
00124 // pixel lies directly between two components and cannot be immediately 
00125 // assigned.  A masked pixel is not inside the targeted region of the
00126 // sub-componentmap and is not used in decomposition or fitting.
00127    enum componentValues {
00128      INDETERMINATE = -1,
00129      MASKED = -2 
00130    };
00131 
00132    ImageDecomposer() = delete;
00133 
00134    ImageDecomposer(const ImageInterface<T>& image);
00135 
00136    ImageDecomposer(const ImageDecomposer<T>& other);
00137 
00138    ImageDecomposer<T> &operator=(const ImageDecomposer<T> &other) = delete;
00139 
00140    ~ImageDecomposer();
00141 
00142 // Tell the decomposer what image to decompose ("target image").
00143 // Also resets the internal component map.
00144 //   void setImage (ImageInterface<T>& image);
00145 
00146 // Tells the program whether or not to use the contour-based deblender. If not,
00147 // the program will instead perform a single thresholding followed by a
00148 // local maximum scan before fitting.
00149    void setDeblend(Bool deblendIt=True);
00150 
00151 // Specifies deblending options:
00152 // <ul>
00153 // <li> thresholdVal: noise cutoff level, used to distinguish source pixels
00154 //      from background pixels.  Also, regions which are not blended above this
00155 //      value will be fit separately.
00156 // <li> nContour: number of total contours to use in deblending regions.
00157 // <li> minRange: the minimum number of contours necessary to distinguish
00158 //      an object as a separate component.
00159 // <li> nAxis: paramater used to define whether or not two adjacent blocks
00160 //      of pixels are contiguous - see identifyRegions for details.
00161 // </ul>
00162 // See decomposeImage for more information on the deblending process.
00163    void setDeblendOptions(T thresholdVal=0.1, uInt nContour=11, 
00164                           Int minRange=2, Int nAxis=2);
00165 
00166 // Tells the program whether or not to perform fitting.  If not, the component
00167 // list will be dstermined by estimation from the values of the first and 
00168 // second order moments of each component.
00169    void setFit(Bool fitIt=True);
00170 
00171 // Specifies fitting options:
00172 // <ul>
00173 // <li> maximumRMS: The maximum RMS residual value (after fitting) allowed to
00174 //      identify a fit as successful.
00175 // <li> maxRetries: the maximum number of times the fit will be restarted in
00176 //      order to try to reach a successful result (convergent with 
00177 //      RMS < maximumRMS).  The default value of -1 tells the program
00178 //      to calculate a reasonable value automatically based on the complexity
00179 //      of the image.
00180 // <li> maxIter: maximum number of iterations to spend on each fit.
00181 // <li> convCriteria: criterion to establish convergence: see NonLinearFitLM.
00182 // </ul>
00183 // Additional information on these parameters can be found in FitGaussian.
00184    void setFitOptions(T maximumRMS=0.1, Int maxRetries=-1, uInt maxIter=256,
00185                       T convCriteria=0.0001);
00186 
00187 
00188 // The primary method of this class - executes the instructions stated in the
00189 // options above by deblending and/or fitting to the image to generate
00190 // the component map and/or component list.
00191   void decomposeImage(); 
00192 
00193 
00194 // Returns the number of regions found in the image.  A 'region' as defined
00195 // in this code is a subset of the image of contiguous pixels whose values
00196 // are greater than the threshold value specified in decomposeImage. A region
00197 // may contain one or more components.
00198   uInt numRegions() const;
00199 
00200 // Returns the number of components found in the image.  A 'component' as
00201 // defined in this code is a source that can be described as a single Gaussian.
00202 // This can only be determined after deblending.
00203   uInt numComponents() const;
00204 
00205 // Returns the shape of the component map.
00206   IPosition shape() const;                
00207 
00208 // Returns the length of a specific axis.
00209   Int shape(uInt axis) const;            
00210 
00211 // Returns True if the image has been thresholded (split up into regions.)
00212   Bool isDerived() const;
00213 
00214 // Returns True if the image has been decomposed (split up into components.)
00215   Bool isDecomposed() const;  
00216 
00217 // Returns the component parameters as a Matrix.  (Ideally, this should be
00218 // a ComponentList.)
00219   Matrix<T> componentList() const;
00220 
00221 // Currently does nothing; in the future should return the component map
00222 // in a way that it can be seen by the user in AIPS++, preferably as a
00223 // colorized image.
00224   void componentMap() const; 
00225 
00226 // Command-line text output functions.
00227 // <group>
00228   void display() const;
00229   void displayContourMap(const Vector<T>& clevels) const;
00230   void printComponents() const;
00231 // </group>
00232 // Boxes each region in the componentmap:
00233 // blc is set to the lowest coordinate value in each region;
00234 // trc is set to one above the highest coordinate value in each region.
00235   void boundRegions(Block<IPosition>& blc, Block<IPosition>& trc);
00236 private:
00237   ImageInterface<T> *itsImagePtr;// Points to the target image.
00238   Lattice<Int> *itsMapPtr;       // The actual component map.  
00239   IPosition itsShape;            // Component map shape
00240   uInt itsDim;                   // Component map number of dimensions
00241   uInt itsNRegions;              // Number of distinct regions in component map
00242   uInt itsNComponents;           // Number of components that have been fitted
00243   Matrix<T> itsList;             // The component list (Gaussian parameters for
00244                                  // each component.)  
00245   Bool itsDeblendIt;
00246   T itsThresholdVal;
00247   uInt itsNContour;              // Decomposition options
00248   Int itsMinRange;               // IMPR: maybe use a struct?
00249   Int itsNAxis;
00250   
00251   Bool itsFitIt;
00252   T itsMaximumRMS; 
00253   Int itsMaxRetries;             // Fitting options
00254   uInt itsMaxIter;
00255   T itsConvCriteria;
00256 
00257 
00258 
00259   void copyOptions(const ImageDecomposer<T>& other);
00260 
00261 // Makes sure a pair of IPositions is in the correct format for blc/trc, and
00262 // corrects them if they are not.
00263   void correctBlcTrc(IPosition& blc, IPosition& trc) const;
00264 
00265 // Used as an N-dimensional interator.  This should probably be replaced by
00266 // LatticeIterators...?
00267 // <group>
00268   Bool increment(IPosition& pos, const IPosition& shape) const;
00269   void decrement(IPosition& pos) const;
00270 // </group>
00271 
00272 // Returns the component to which the specified cell belongs
00273 // <group>
00274   Int getCell(Int x, Int y) const;             
00275   Int getCell(Int x, Int y, Int z) const;    
00276   Int getCell(const IPosition& coord) const;
00277 // </group>
00278 
00279 // Assigns the specified cell to the specified component
00280 // <group>
00281   void setCell(Int x, Int y, Int sval);       
00282   void setCell(Int x, Int y, Int z, Int sval); 
00283   void setCell(const IPosition& coord, Int sval);
00284 // </group>
00285 
00286 // Semi-automatic way to set contour levels: at the given increment counting
00287 // between mincon and maxcon.
00288   Vector<T> autoContour(T minCon, T maxCon, T inc) const;
00289 
00290 // Linearly spaces contours between minvalue and just below the
00291 // maximum value in the target region of the target image, and returns
00292 // the contour values as a Vector.
00293   Vector<T> autoContour(Int nContours=11, T minValue=0) const;
00294 
00295 // Nonlinear spacing option for contouring; spaces contours according to the
00296 // function given.  The domain of the function is 0 <-> ncontours-1; the
00297 // range is automatically calibrated to be minvalue <-> maxvalue.  The function
00298 // should be nondecreasing in the domain such that each contour is greater
00299 // than the last.
00300   Vector<T> autoContour(const Function1D<T>& fn,
00301                         Int nContours = 11, T minValue = 0) const;
00302 
00303 
00304 //Eliminates any regions whose corresponding values in killRegion are True
00305 // by setting all pixel values in the componentmap set to that region to
00306 // zero.  Zero-oriented; there is an offset of one between the index in
00307 // killRegion and the actual region in the componentmap.
00308   void destroyRegions(const Vector<Bool>& killRegion);
00309 
00310 // Eliminates regions with no cells by replacing them with higher-numbered 
00311 // regions.
00312   void renumberRegions();
00313 
00314 // Overlays a smaller map onto an empty region of a larger map,
00315 // and adds submap component list to main component list.
00316 // The user should exercise caution with this function and synthesize submaps
00317 // only into regions of the main map that are truly empty (0), as no blending
00318 // is assumed between different maps.
00319   void synthesize(const ImageDecomposer<T>& subdecomposer, IPosition blc);
00320 
00321 // Set all elements in the component map to zero and clear the component list.
00322   void zero();
00323 
00324 // Set all nonmasked elements in the component map to zero and clear the 
00325 // component list.
00326   void clear();
00327 
00328 
00329 
00330 // Finds the greatest value inside the specified rectangular area of the
00331 // target image.
00332 // <group>
00333   T findAreaGlobalMax(IPosition blc, IPosition trc) const;
00334   void findAreaGlobalMax(T& maxval, IPosition& maxvalpos, 
00335                          IPosition blc, IPosition trc) const;
00336   Vector<T> findAreaGlobalMax(IPosition blc, IPosition trc, Int naxis) const;
00337   void findAreaGlobalMax(Vector<T>& maxvals,
00338                          Block<IPosition>& maxvalpos,
00339                          IPosition blc, IPosition trc,
00340                          Int naxis) const;
00341 // </group>
00342 
00343 // Finds all local maxima inside the specified rectangular area of the
00344 // target image.
00345 // <group>
00346   Vector<T> findAreaLocalMax(IPosition blc, IPosition trc, Int naxis) const;
00347   void findAreaLocalMax(Vector<T>& maxvals,Block<IPosition>& maxvalpos, 
00348                         IPosition blc, IPosition trc, Int naxis) const;
00349 // </group>
00350 
00351 // Finds the maximum value of the target image in each region of the
00352 // componentmap.
00353 // <group>
00354   Vector<T> findAllRegionGlobalMax() const;
00355   void findAllRegionGlobalMax(Vector<T>& maxvals, 
00356                               Block<IPosition>& maxvalpos) const;
00357 // </group>
00358 
00359 
00360 // Finds all local maxima of the target image inside the specifed region
00361 // of the componentmap.
00362 // <group>
00363   Vector<T> findRegionLocalMax(Int nregion, Int naxis) const;
00364   void findRegionLocalMax(Vector<T>& maxvals, Block<IPosition>& maxvalpos, 
00365                           Int nregion, Int naxis) const;
00366 // </group>
00367 
00368 
00369 //Compares specified pixel to adjacent pixels to determine if it is
00370 //greatest in local pixel block.
00371 //2D:
00372 //naxis = 1: compare to 4 adjacent pixels (axes only)
00373 //naxis = 2: compare to 8 adjacent pixels (axes and diagonals)
00374 //3D:
00375 //naxis = 1: compare to 6 adjacent pixels (axes only)
00376 //naxis = 2: compare to 18 adjacent pixels (axes and 2-axis diagonals)
00377 //naxis = 3: compare to 26 adjacent pixels (axes and 2/3-axis diagonals)
00378 // <group>
00379   Bool isLocalMax(const IPosition& pos, Int naxis) const;
00380   Bool isLocalMax(Int x, Int y, Int naxis) const;
00381   Bool isLocalMax(Int x, Int y, Int z, Int naxis) const;
00382 // </group>
00383 
00384 // Finds a rough estimate of the width of each component by scanning to find
00385 // the full width at quarter maximum.
00386 // Requires the location of each component.
00387 // This function is mostly obsolete, and is only used when the contour 
00388 // deblender is off (since the component map is necessary to determine the
00389 // moments).
00390   void estimateComponentWidths(Matrix<T>& width,
00391                                const Block<IPosition>& maxvalpos) const;
00392 
00393 // Calculates the 0th-2nd order moments of a region.
00394   Array<T> calculateMoments(Int region) const;
00395 
00396 
00397 // Performs a single threshold scan on the image.  In other words,
00398 // identifies all contigous blocks of pixels in the target image above the
00399 // threshold value thrval, assigning each unique block to an integer,
00400 // starting at one.  All pixels with target image values below thrval are set
00401 // to zero.
00402   uInt identifyRegions(T thrval, Int naxis=2);
00403 
00404 // Performs the contour decomposition on a blended image to generate a 
00405 // component map that can detect components blended above any threshold(s),
00406 // by performing threshold scans at each contour level and recognizing
00407 // as individual any components that are distinct above any such level.
00408   void deblendRegions(const Vector<T>& contours, Int minRange=1, Int naxis=2);
00409 
00410 // Retrieves the target image's value at the given location.
00411 // <group>
00412   T getImageVal(IPosition coord) const;
00413   T getImageVal(Int x, Int y) const;
00414   T getImageVal(Int x, Int y, Int z) const;
00415 // </group>
00416 
00417 // Retrieves the number of the highest contour with a value less then the
00418 // target image's value at the given location.
00419 // <group>
00420   Int getContourVal(IPosition coord, const Vector<T>& clevels) const;
00421   Int getContourVal(Int x, Int y, Int z, const Vector<T>& clevels) const;
00422   Int getContourVal(Int x, Int y, const Vector<T>& clevels) const;
00423   Int getContourVal(T val, const Vector<T>& clevels) const;
00424 // </group>
00425 
00426 // Fits multiple gaussians to a single region.  First performs a local 
00427 // maximum scan to estimate the number of components in the region.
00428   Matrix<T> fitRegion(Int region);
00429 
00430 // Fits gaussians to an image; multiple gaussians per region in the component
00431 // map.  The regions are fit sequentially and independently, so this function 
00432 // can be used on the main image.  If the map is not yet thresholded, will fit
00433 // to the entire image as if it  were a single composite object, which will be
00434 // very slow.
00435   void fitRegions();
00436 
00437 // Fits gaussians to an image; one gaussian per region in the pmap.
00438 // This function is intended to be used only by ImageDecomposer on its
00439 // intermediary subimages; using it at higher level will execute a full
00440 // gaussian fit on the main image and will be extremely slow. Every 
00441 // nonflagged object pixel in the image is used in fitting.
00442 
00443 // If the deblended flag is True, the function will treat each region as
00444 // an individual component and will fit that many gaussians to the image
00445   void fitComponents();
00446 
00447 // Estimate the component parameters based on moments calculated using 
00448 // the component map.
00449   Matrix<T> estimateComponents();
00450 
00451 // Fits the specified number of 3D gaussians to the data, and returns 
00452 // solution in image (world) coordinates.  Essentially just an interface
00453 // for FitGaussian.
00454 // <group>
00455   Matrix<T> fitGauss(const Matrix<T>& positions, const Vector<T>& dataValues,
00456                      const Matrix<T>& initestimate) const;
00457 
00458 // </group>
00459 
00460 };
00461 
00462 } //# NAMESPACE CASA - END
00463 
00464 #ifndef CASACORE_NO_AUTO_TEMPLATES
00465 #include <imageanalysis/ImageAnalysis/ImageDecomposer.tcc>
00466 #endif //# CASACORE_NO_AUTO_TEMPLATES
00467 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

Generated on 31 Aug 2016 for casa by  doxygen 1.6.1