ImageDecomposer.h

ImageDecomposer is an image decomposition tool that performs several tasks, with the end result being that a strongly blended image is separated into components - both in the sense that it determines the parameters for each component (assuming a Gaussian model) and that it physically assigns each pixel in the image to an individual object. The products of these two operations are called the component list and the component map, respectively. The fitting process (which determines the component list) and the pixel-decomposition process (which determines the component map) are designed to work cooperatively to increase the efficiency and accuracy of both, though each can operate without the other if necessary.

The algorithm between the decomposition is based on the function clfind described in Williams et al 1994, which uses a contouring procedure whereby a closed contour designates a separate component. The program first separates the image into clearly distint 'regions' of blended emission, then contours each region to determine the areas constituting each component and passes this information on to the fitter, which determines the component list.

The software is compatible with 2 and 3 dimensional images, but is not yet structured for higher dimensions.

Example

  TempImage<Double> image;
  //(populate the image with data: see dImageDecomposer.cc)
  ImageDecomposer<Double> id(image);
  id.setDeblendOptions(0.3, 8);
  id.setFitOptions(0.4);
  id.decomposeImage();
  id.display();
  id.printComponents();

Motivation

To Do

Generalize dimensionality
Use Lattice iterators in place of IPosition loops wherever possible
Speed up fitting by not sending every region pixel to the fitter
Send the completed componentmap to the user as an AIPS++ (Int?) Image
Return a ComponentList instead of a Matrix
Enable custom contouring at user level
Add progress meter
Numerous other improvements to make are documented in the code

Member Description

enum componentValues

'Special' flag values for pixels in the component map. An indeterminate pixel lies directly between two components and cannot be immediately assigned. A masked pixel is not inside the targeted region of the sub-componentmap and is not used in decomposition or fitting.

ImageDecomposer()

Default constructor. Object is not viable until setImage called

ImageDecomposer(ImageInterface<T>& image)

Construct from image

ImageDecomposer(const ImageDecomposer<T>& other)

Copy constructor.

ImageDecomposer<T> &operator=(const ImageDecomposer<T> &other)

Assignment

~ImageDecomposer()

Destructor

void setImage (ImageInterface<T>& image)

Tell the decomposer what image to decompose ("target image"). Also resets the internal component map.

void setDeblend(Bool deblendIt=True)

Tells the program whether or not to use the contour-based deblender. If not, the program will instead perform a single thresholding followed by a local maximum scan before fitting.

void setDeblendOptions(T thresholdVal=0.1, uInt nContour=11, Int minRange=2, Int nAxis=2)

Specifies deblending options:

thresholdVal: noise cutoff level, used to distinguish source pixels from background pixels. Also, regions which are not blended above this value will be fit separately.
nContour: number of total contours to use in deblending regions.
minRange: the minimum number of contours necessary to distinguish an object as a separate component.
nAxis: paramater used to define whether or not two adjacent blocks of pixels are contiguous - see identifyRegions for details.

See decomposeImage for more information on the deblending process.

void setFit(Bool fitIt=True)

Tells the program whether or not to perform fitting. If not, the component list will be dstermined by estimation from the values of the first and second order moments of each component.

void setFitOptions(T maximumRMS=0.1, Int maxRetries=-1, uInt maxIter=256, T convCriteria=0.0001)

Specifies fitting options:

maximumRMS: The maximum RMS residual value (after fitting) allowed to identify a fit as successful.
maxRetries: the maximum number of times the fit will be restarted in order to try to reach a successful result (convergent with RMS < maximumRMS). The default value of -1 tells the program to calculate a reasonable value automatically based on the complexity of the image.
maxIter: maximum number of iterations to spend on each fit.
convCriteria: criterion to establish convergence: see NonLinearFitLM.

Additional information on these parameters can be found in FitGaussian.

void decomposeImage()

The primary method of this class - executes the instructions stated in the options above by deblending and/or fitting to the image to generate the component map and/or component list.

uInt numRegions() const

Returns the number of regions found in the image. A 'region' as defined in this code is a subset of the image of contiguous pixels whose values are greater than the threshold value specified in decomposeImage. A region may contain one or more components.

uInt numComponents() const

Returns the number of components found in the image. A 'component' as defined in this code is a source that can be described as a single Gaussian. This can only be determined after deblending.

IPosition shape() const

Returns the shape of the component map.

Int shape(uInt axis) const

Returns the length of a specific axis.

Bool isDerived() const

Returns True if the image has been thresholded (split up into regions.)

Bool isDecomposed() const

Returns True if the image has been decomposed (split up into components.)

Matrix<T> componentList() const

Returns the component parameters as a Matrix. (Ideally, this should be a ComponentList.)

void componentMap() const

Currently does nothing; in the future should return the component map in a way that it can be seen by the user in AIPS++, preferably as a colorized image.

void display() const
void displayContourMap(const Vector<T>& clevels) const
void printComponents() const

Command-line text output functions.

void copyOptions(const ImageDecomposer<T>& other)

void correctBlcTrc(IPosition& blc, IPosition& trc) const

Makes sure a pair of IPositions is in the correct format for blc/trc, and corrects them if they are not.

Bool increment(IPosition& pos, const IPosition& shape) const
void decrement(IPosition& pos) const

Used as an N-dimensional interator. This should probably be replaced by LatticeIterators...?

Int getCell(Int x, Int y) const
Int getCell(Int x, Int y, Int z) const
Int getCell(const IPosition& coord) const

Returns the component to which the specified cell belongs

void setCell(Int x, Int y, Int sval)
void setCell(Int x, Int y, Int z, Int sval)
void setCell(const IPosition& coord, Int sval)

Assigns the specified cell to the specified component

Vector<T> autoContour(T minCon, T maxCon, T inc) const

Semi-automatic way to set contour levels: at the given increment counting between mincon and maxcon.

Vector<T> autoContour(Int nContours=11, T minValue=0) const

Linearly spaces contours between minvalue and just below the maximum value in the target region of the target image, and returns the contour values as a Vector.

Vector<T> autoContour(const Function1D<T>& fn, Int nContours = 11, T minValue = 0) const

Nonlinear spacing option for contouring; spaces contours according to the function given. The domain of the function is 0 <-> ncontours-1; the range is automatically calibrated to be minvalue <-> maxvalue. The function should be nondecreasing in the domain such that each contour is greater than the last.

void ImageDecomposer<T>::destroyRegions(const Vector<Bool>& killRegion)

Eliminates any regions whose corresponding values in killRegion are True by setting all pixel values in the componentmap set to that region to zero. Zero-oriented; there is an offset of one between the index in killRegion and the actual region in the componentmap.

void renumberRegions()

Eliminates regions with no cells by replacing them with higher-numbered regions.

void synthesize(const ImageDecomposer<T>& subdecomposer, IPosition blc)

Overlays a smaller map onto an empty region of a larger map, and adds submap component list to main component list. The user should exercise caution with this function and synthesize submaps only into regions of the main map that are truly empty (0), as no blending is assumed between different maps.

void zero()

Set all elements in the component map to zero and clear the component list.

void clear()

Set all nonmasked elements in the component map to zero and clear the component list.

void boundRegions(Block<IPosition>& blc, Block<IPosition>& trc)

Boxes each region in the componentmap: blc is set to the lowest coordinate value in each region; trc is set to one above the highest coordinate value in each region.

T findAreaGlobalMax(IPosition blc, IPosition trc) const
void findAreaGlobalMax(T& maxval, IPosition& maxvalpos, IPosition blc, IPosition trc) const
Vector<T> findAreaGlobalMax(IPosition blc, IPosition trc, Int naxis) const
void findAreaGlobalMax(Vector<T>& maxvals, Block<Block>& maxvalpos, Block blc, Block trc, Int naxis) const

Finds the greatest value inside the specified rectangular area of the target image.

Vector<T> findAreaLocalMax(IPosition blc, IPosition trc, Int naxis) const
void findAreaLocalMax(Vector<T>& maxvals,Block<Block>& maxvalpos, Block blc, Block trc, Int naxis) const

Finds all local maxima inside the specified rectangular area of the target image.

Vector<T> findAllRegionGlobalMax() const
void findAllRegionGlobalMax(Vector<T>& maxvals, Block<Block>& maxvalpos) const

Finds the maximum value of the target image in each region of the componentmap.

Vector<T> findRegionLocalMax(Int nregion, Int naxis) const
void findRegionLocalMax(Vector<T>& maxvals, Block<Block>& maxvalpos, Int nregion, Int naxis) const

Finds all local maxima of the target image inside the specifed region of the componentmap.

Bool isLocalMax(const IPosition& pos, Int naxis) const
Bool isLocalMax(Int x, Int y, Int naxis) const
Bool isLocalMax(Int x, Int y, Int z, Int naxis) const

Compares specified pixel to adjacent pixels to determine if it is greatest in local pixel block. 2D: naxis = 1: compare to 4 adjacent pixels (axes only) naxis = 2: compare to 8 adjacent pixels (axes and diagonals) 3D: naxis = 1: compare to 6 adjacent pixels (axes only) naxis = 2: compare to 18 adjacent pixels (axes and 2-axis diagonals) naxis = 3: compare to 26 adjacent pixels (axes and 2/3-axis diagonals)

void estimateComponentWidths(Matrix<T>& width, const Block<Block>& maxvalpos) const

Finds a rough estimate of the width of each component by scanning to find the full width at quarter maximum. Requires the location of each component. This function is mostly obsolete, and is only used when the contour deblender is off (since the component map is necessary to determine the moments).

Array<T> calculateMoments(Int region) const

Calculates the 0th-2nd order moments of a region.

uInt identifyRegions(T thrval, Int naxis=2)

Performs a single threshold scan on the image. In other words, identifies all contigous blocks of pixels in the target image above the threshold value thrval, assigning each unique block to an integer, starting at one. All pixels with target image values below thrval are set to zero.

void deblendRegions(const Vector<T>& contours, Int minRange=1, Int naxis=2)

Performs the contour decomposition on a blended image to generate a component map that can detect components blended above any threshold(s), by performing threshold scans at each contour level and recognizing as individual any components that are distinct above any such level.

T getImageVal(IPosition coord) const
T getImageVal(Int x, Int y) const
T getImageVal(Int x, Int y, Int z) const

Retrieves the target image's value at the given location.

Int getContourVal(IPosition coord, const Vector<T>& clevels) const
Int getContourVal(Int x, Int y, Int z, const Vector<T>& clevels) const
Int getContourVal(Int x, Int y, const Vector<T>& clevels) const
Int getContourVal(T val, const Vector<T>& clevels) const

Retrieves the number of the highest contour with a value less then the target image's value at the given location.

Matrix<T> fitRegion(Int region)

Fits multiple gaussians to a single region. First performs a local maximum scan to estimate the number of components in the region.

void fitRegions()

Fits gaussians to an image; multiple gaussians per region in the component map. The regions are fit sequentially and independently, so this function can be used on the main image. If the map is not yet thresholded, will fit to the entire image as if it were a single composite object, which will be very slow.

void fitComponents()

Fits gaussians to an image; one gaussian per region in the pmap. This function is intended to be used only by ImageDecomposer on its intermediary subimages; using it at higher level will execute a full gaussian fit on the main image and will be extremely slow. Every nonflagged object pixel in the image is used in fitting.

If the deblended flag is True, the function will treat each region as an individual component and will fit that many gaussians to the image

Matrix<T> estimateComponents()

Estimate the component parameters based on moments calculated using the component map.

Matrix<T> fitGauss(const Matrix<T>& positions, const Vector<T>& dataValues, const Matrix<T>& initestimate) const

Fits the specified number of 3D gaussians to the data, and returns solution in image (world) coordinates. Essentially just an interface for FitGaussian.