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