DynamicPSF

PSF Model Functions

The current version of DynamicPSF can fit elliptical Gaussian and Moffat ^[1] functions. These functions have been selected because their shapes make them particularly suitable to model stellar objects represented on most deep-sky images.

DynamicPSF implements the Levenberg-Marquadt algorithm as a model trust region minimization procedure. Given sampled data consisting of two sets of independent and dependent variables, {x₀ … x_N–1} and {y₀ … y_N–1} respectively, the algorithm minimizes the sum of absolute differences of the model function f with the sampled data for a vector a of fitted function parameters:

In this application the independent variables are pairs of pixel coordinates on the image, and the dependent variables are the corresponding pixel values. The underlying implementation is a custom adaptation of the MINPACK library ^{[2] [3]} to the PixInsight/PCL framework. Some routines also implement ideas and algorithms from other sources listed in the References section.^{[4] [5]}

Gaussian Function

An elliptical Gaussian PSF fitting function is defined as:

where the parameters are as follows:

• B — Average local background.
• A — Amplitude, which is the maximum value of the fitted PSF, and also the function's value at the centroid coordinates.
• x_0, y₀ — Centroid coordinates in pixel units. This is the position of the center of symmetry of the fitted PSF.
• σ_x, σ_y — Standard deviation of the Gaussian distribution on the horizontal and vertical axes, measured in pixels.

Moffat Function

An elliptical Moffat PSF fitting function is defined as:

where the parameters are as in the Gaussian case, and the β exponent controls the overall shape of the fitting function. In the current implementation of DynamicPSF, the β parameter can vary in the range ]0,10]. When β=1, the above equation corresponds to a Lorentzian distribution.

Rotated Functions

DynamicPSF also supports rotated PSF fitting functions. When the difference between σ_x and σ_y is larger than or equal to 0.01 pixel (the nominal fitting resolution), DynamicPSF fits an additional θ parameter, which is the rotation angle of the X axis with respect to the centroid coordinates. θ varies in the range [0°,180°[. For a rotated PSF, the x and y coordinates in the equations above must be replaced by their rotated counterparts x' and y', respectively:

Rotation angles are measured in counter-clockwise direction, and can be represented either as unsigned values in the [0°,180°[ range, or as signed quantities in the [–90°,+90°] range. In all cases the functions are fitted to ensure that σ_x ≥ σ_y.

Graphical Function Representations

Graphical representations are the best way to understand the geometries of different PSF fitting functions and how they change by varying their parameters. The figure below represents several functions graphically. Two-dimensional representations have been generated by running the following script in PixInsight:

#include <pjsr/UndoFlag.jsh>

function filterToImage( filter, id )
{
   var w = new ImageWindow( 1, 1 );
   var v = w.mainView;
   v.beginProcess( UndoFlag_NoSwapFile );
   v.image.assign( filter.toImage() );
   v.id = id;
   v.endProcess();
   w.show();
   w.zoomToFit();
}

#define FILTER_SIZE  101

filterToImage( Matrix.gaussianFilterBySize( FILTER_SIZE ), "Gaussian" );
filterToImage( Matrix.moffatFilterBySize( FILTER_SIZE, 1.5 ), "Moffat_15" );
filterToImage( Matrix.moffatFilterBySize( FILTER_SIZE, 4.0 ), "Moffat_40" );
filterToImage( Matrix.gaussianFilterBySize( FILTER_SIZE, 0.01,
                                            0.5, Math.rad( 45 ) ), "GaussianRotated" );

The script represents four PSF fitting functions (those included in the following figure) as new image windows. This script can be easily modified to represent other Gaussian and Moffat functions with different parameters. Once we have the functions represented as images, the standard 3DPlot script can be used in PixInsight to generate high-quality, three-dimensional renditions.

Figure 1

Gaussian function.

Moffat function, β=1.5.

Moffat function, β=4.

Gaussian function, σ_x = 2σ_y, θ = 45°.

Full Width at Half Maximum (FWHM)

For all PSF model functions DynamicPSF computes two additional full width at half maximum (FWHM) values:

• FWHM_x — The FWHM on the X axis.
• FWHM_y — The FWHM on the Y axis. For circular (undistorted, unrotated) functions, we have FWHM_y = FWHM_x.

The FWHM is a well-known, standardized and easy to understand measurement of the size of a star as seen on the image. It is the width in pixels of the fitted function, measured horizontally at half its height. The FWHM on the X axis is given by:

respectively for Gaussian and Moffat fitting functions, with similar expressions for the Y axis by replacing σ_x with σ_y. Note that FWHM measurements for Gaussian and Moffat functions are in general not compatible and should not be mutually compared.

Goodness of Fit

For each fitted PSF DynamicPSF provides an estimate of the goodness of fit, or how well the computed PSF function agrees with the actual image pixels that have been sampled to perform the function fitting process. In the current version of DynamicPSF a mean absolute difference (MAD) function is computed for this purpose:

where f represents the fitted function and the I_i elements are the sampled pixels from the original image. The MAD function is intended to provide a robust estimate of the suitability of a model function and its fitted parameters to represent the actual sampled data. MAD estimates should be used to remove outliers from the set of fitted objects during a DynamicPSF session. In this context, outliers correspond to poorly sampled objects, such as saturated or too dim objects, or objects of the wrong type (e.g. nonstellar objects such as small background galaxies) that can contaminate the computed average PSF.

PSF Data

The first and largest control on the DynamicPSF interface is a tree box control where the data for all fitted PSF model functions are shown in a hierarchical way. In the first hierarchy level we have PSF collections. There is one PSF collection item for each image involved in the current DynamicPSF session. Each PSF collection can have an unlimited number of star item descendants. Each star item, in turn, has one or more PSF fits performed for the same object, whose parameters are listed in the corresponding rows of the tree box. For color images, each channel is treated as a separate image, so there can be one, two or three PSF fits performed with the same parameters for each nominal channel of a RGB color image.

The following table describes all the data items provided by the current version of DynamicPSF.

For circular functions, the sy, FWHMy, r and theta data items are not provided. For Gaussian functions the beta item is not provided.

Creating New PSF Fits

A new PSF fitting process is performed when you click on an image and release the pointing device without dragging. If there are no PSF fits for the image where you clicked, a new PSF collection is created. A new star item is then created and appended at the end of the corresponding PSF collection, either a newly created or an already existing collection. Note that, depending on the selected options in PSF Model Functions section, this operation can create one or several PSF fits for the same object. As we have said before, you can repeat this operation unlimited times.

When you click on an image to create a new PSF fit, DynamicPSF automatically looks for the nearest small-scale, bright and concentrated image structure, so you don't have to click at any exact location; just clicking around a star is sufficient. The object detection routine is both robust and stable. It doesn't matter where you click on the image; as far as you click reasonably close to an object of interest, the routine will always converge at the same pixel coordinates and the fitted PSF parameters will always be the same for a given object, irrespective of the initial search location.

Figure 2

For each PSF fit, DynamicPSF generates a set of vector graphics elements on the target image. The PSF coordinate system is rotated by the angle θ with respect to the image coordinate system. In a three-dimensional representation of the PSF model function, the Z axis is perpendicular to the image plane (the XY plane), and Z coordinates correspond to image pixel values. The plotted ellipse represents the intersection between the 3D PSF model function and a plane parallel to the image plane, elevated by half the maximum value of the function. This defines the full width at half maximum on the X and Y directions, FWHMx and FWHMy, respectively.

Working with Color images

When working with color images, three PSF fits are created by default for each fitted star, one for each nominal RGB channel (alpha channels are always ignored). You can restrict creation of new PSF fits for a particular channel by selecting it for display on the target image. For example, if you select the green channel (using the channel selector combo box, or pressing Ctrl+Shift+G), only one PSF fit will be created for that channel.

Selecting PSF Fits

You can select existing PSF fits in several ways:

• By clicking within an existing PSF fit on the image. Each PSF fit is plotted as an ellipse inside a bounding rectangle. To select a PSF fit you have to click inside its bounding rectangle. Hold the Shift key pressed while you click to add a PSF fit to the current selection.
• By clicking on the image and dragging to define a selection rectangle. All existing PSF fits within the selection rectangle will be selected.
• Using the tree box control on the DynamicPSF interface. You can click on a row of the tree box to select the corresponding star item and/or associated PSF fit(s), Shift+click or click and drag to define multiple selections, etc.

Selected PSF fits are plotted with a special color on the image (green by default). Note that the current selection is always coherent between the vector graphics plotted on the image(s) and the items in the DynamicPSF interface.

Removing Items

PSF fits can be removed in two ways:

• By clicking the Remove button. This will delete all selected objects, which can be star items (each item with one or more PSF fits) and/or individual PSF fits, depending on the current selection.
• By defining a selection and pressing the Del or Backspace key. This works in the same way as the Remove button described above.

Be careful as item removal operations cannot be undone.

Selecting PSF Model Functions

The PSF Model Functions section of the DynamicPSF interface allows you to select one or more PSF model functions to use in new PSF fits. The following options are available:

Auto

Select this option to let DynamicPSF find the best model function for each fitted object. When this option is selected, DynamicPSF fits a Gaussian function and one or several Moffat functions, with variable or fixed β exponent, and selects the function that best reproduces the sampled image data. This is the default option. Note that when this option is selected, the Gaussian and Moffat options (see below) are unavailable.

Gaussian

Select this option to fit Gaussian PSF model functions.

Moffat

Select this option to fit Moffat PSF model functions with variable (that is, fitted) β parameters.

Moffat10, Moffat8, Moffat6, Moffat4, Moffat25, Moffat15, Lorentzian

Select one or more of these options to fit Moffat PSF model functions with fixed β parameters of 10, 8, 6, 4, 2.5, 1.5 and 1.0, respectively (a Lorentzian function is equivalent to a Moffat function with β=1).

Circular PSF

Select this option to fit circular PSF functions. Disable it to fit elliptical functions. An elliptical function provides two separate axes and a rotation angle. In general, elliptical functions are preferable because they provide more information about the true shapes and orientations of the fitted PSFs. Sometimes, however, circular functions may be preferable. For example, strongly undersampled images usually don't provide enough data to fit elliptical functions reliably. Uncertainty due to high noise levels also affects complex function fittings adversely. In such cases a circular function can provide more robust and useful results. Elliptical PSF functions are fitted by default.

Signed angles

When this option is enabled, rotation angles are shown as signed values in the [–90°,+90°;] range. When this option is disabled, rotation angles are represented as unsigned quantities in the [0°,180°[ range. Signed angles are useful to prevent ambiguities introduced by small rotations around zero degrees. These 'oscillatory' rotations occur frequently due to uncertainty in fitting rotation angles for nearly circular stars. For example, imagine a set of PSF fits where we have two subsets with rotations around zero and 180°. If we compute average PSF parameters, the resulting average rotation would be around 90°. This may be incorrect if both subsets are actually due to dispersion caused by fitting uncertainty, in which case the correct average rotation would be a value close to zero degrees. This option is enabled by default.

Star Detection Parameters

The Star Detection section provides options to control how DynamicPSF searches and detects stellar objects for PSF fitting:

Search radius

This parameter determines the size in pixels of the initial search box used to detect stars when you clic on an image. Increase it to favor detection of larger stars. Decrease it to facilitate selection of smaller stars. For example, a smaller search radius may be necessary to deal with dense star fields. The valid range is from one to 127 pixels. The default value of eight pixels is quite appropriate in most cases.

Background threshold

Threshold value for rejection of background pixels, in sigma units. This value is used by the object detection routine to classify sampled pixels into two disjoint sets: pixels belonging to the local background, and pixels belonging to the object being detected. A smaller threshold value makes the star detection routine less sensitive. This can be useful to avoid detection of very faint objects. On the other hand, a higher threshold allows you to isolate very small and faint objects, which may or may not be a good idea, depending on what you want to do and on the image you are working with. This parameter can range from 0.05 to 5 in sigma units. The default value of one sigma is normally appropriate in most cases.

Automatic aperture

A PSF fitting process uses a square sampling region around the fitted object to gather a set of source image pixels. Then a PSF model function is fitted to minimize the difference with the sampled source pixels. When this option is enabled, DynamicPSF computes optimal sampling dimensions adaptively for each fitted star.

If the sampling region is too small, the local background cannot be determined accurately, and hence the PSF fit will be less accurate. If the sampling region is too large, pixels from nearby objects can contaminate the sample, and computation times grow unnecessarily (bear in mind that they are proportional to the square of the sampling size). The automatic aperture feature of DynamicPSF implements a fast and robust algorithm to find the smallest sampling region necessary to minimize uncertainty in the evaluation of the local background. This option is enabled by default. Disabling it is not recommended under normal working conditions, and can lead to a significant degradation in performance and fitting accuracy.

Image Scale Parameters

This section allows you to specify the scale of the image in arcseconds per pixel. When the image scale is available, DynamicPSF represent FWHM values directly in arcseconds, instead of pixels. The following parameters can be used for this purpose:

Scale mode

Defines a method to compute the image scale. Can be one of:

• Standard FITS keywords. If this mode is selected, DynamicPSF will try to read standard FOCALLEN, XPIXSZ and YPIXSZ FITS header keywords to compute the image scale in arcseconds per pixel. If these keywords are not present or contain invalid values, no image scale is computed and FWHM is expressed in pixels.
• Literal value. This mode allows you to enter the image scale directly in arcseconds per pixel (see the next parameter).
• Custom FITS keyword. This mode allows you to specify the name of a custom FITS keyword whose value is the image scale in arcseconds per pixel.
• Pixels. This will ignore image scale to show FWHM in pixel units. This is the default option.

Image scale

When the literal value mode is selected, this control allows you to enter the image scale in arcseconds per pixel. Note that this value will be used for all the images involved in the current DynamicPSF session.

Custom keyword

When the custom FITS keyword mode is selected, this control allows you to enter the name of the special FITS header keyword that provides the image scale in arcseconds per pixel. If the specified keyword is not available, image scale will be ignored and FWHM will be shown in pixel units.

Tracking Stars

When the Track Stars option is enabled (button pushed), DynamicPSF automatically centers the selected star on its corresponding target image and brings the corresponding view to the top of PixInsight's workspace, if necessary. This is very useful to work at high magnifications, as it frees you from having to zoom in/out and pan the images when you are reviewing a set of PSF fits. This option is disabled by default.

Regenerating PSF Fits

You can regenerate existing PSF fits, which allows you to calculate new PSF parameters after changing some parameters in the PSF Model Functions section. For example, suppose you have generated a number of PSF fits using elliptical functions. Later you decide to fit circular functions for the same stars. You don't need to remove the previous fits and perform the tedious (and error prone) task of redefining new ones; just select the corresponding option (circular functions) and regenerate the PSF fits that you have already defined. You can regenerate in two ways:

• Regenerate selected stars. Click this button to recalculate only the set of selected PSF fits.
• Regenerate all. Click this button to recalculate all PSF fits, irrespective of the current selection.

Sorting PSF Collections

The Sort Stars button allows you to sort all the existing PSF fits in ascending order by selectable criteria. When you click this button, the Sort PSF Collections dialog opens as shown below.

On this dialog you can select one of the following sorting criteria:

• Star Id. Each PSF fit is assigned an identifier, which is a running number that is unique within its parent PSF collection. Sorting by identifier classifies the PSF fits by creation order.
• Background. Sorts PSF fits by the mean local background parameter (B).
• Amplitude. Sorts PSF fits by the amplitude parameter (A), which is the peak value of the fitted PSF model function.
• Standard deviation. Sorts PSF fits by their sizes. Note that the size in the X axis (PSF coordinates) is always used, that is the sx data item. Sorting by this criterion is equivalent to sorting by FWHM.
• Aspect ratio. Sorts by the quotient sy/sx (r data item), that is by the degree of distortion of the fitted PSF model function.
• Rotation angle. Sorts by the rotation angles of fitted PSF model functions (theta parameter). Note that the resulting order depends on whether signed angles are being used in the range [–90°,+90°] or unsigned angles in [0°,180°[. Circular functions are considered to have a zero rotation angle, and are sorted accordingly.
• Rotation angle (absolute value). Same as the preceding criterion, but angles are evaluated in the [0°,180°[ range, irrespective of the current state of the signed angles option.
• Beta exponent. Sorts by the Moffat beta parameter. Gaussian functions are assumed to have a zero beta exponent, and hence are always placed at the top of the list when this sorting criterion is used.
• Mean absolute difference. Sorts by the MAD data item. This is the default sorting criterion. Use it to classify PSF fits by their quality estimates.

Sorting your PSF collections is very useful to find and remove outliers in the sets of fitted stars. The most obvious use of this feature is to sort by fitting quality (the default MAD sorting criterion). MAD estimates allow you to easily cut the lists of PSF fits below a given quality limit. However, other sorting criteria can provide more sophisticated ways to detect incorrect PSF fits that may be difficult to evidence otherwise. For example, if you know that your tracking is accurate and hence all of the stars in your image should have nearly circular shapes, then if there are some objects fitted with relatively small aspect ratios—even if they show pretty low MAD values—, they are clearly nonstellar objects or spurious image features that should be removed from your PSF collections.

Exporting a Synthetic PSF

By clicking the Export Synthetic PSF button, DynamicPSF will generate a new image window with an average PSF computed for the current set of selected PSF fits (you have to select at least one fit for this function to work). Note that the generated PSF image is not the result of averaging all PSF function parameters, but something much more robust and accurate. DynamicPSF generates a synthetic PSF by rendering all fitted PSF model functions superposed by transparency over a black canvas, large enough to accomodate the largest PSF, and all of them centered at their centroid positions. This effectively works as if you had acquired all of the fitted stars centered in the same field of view with multiple exposures of a hypothetical CCD camera, able to accumulate the light from all of them without saturation. The main advantage of this method is that the generated PSF image is almost immune to small variations, and hence very accurate to represent the actual PSF of the image; this quality is what we denote as robustness in this context.

A synthetic PSF is the recommended way to generate a PSF image that can be used with PixInsight's Deconvolution tool as an external PSF. In general, this method provides the best results to deconvolve an image because the obtained PSF is usually very accurate. Achieving similar accuracy by trial and error work with parametric PSFs is a difficult and time consuming task.

Average PSF Parameters

DynamicPSF also allows you to read average PSF parameters for the set of selected PSF fits. This is only possible when the selected PSF fits are congruent, that is when the selected set only contains either Gaussian or Moffat modelling functions. When you click the Average PSF Parameters button the Average Star Data dialog opens and shows a list of parameter values calculated by averaging the set of selected fits, as shown in the following screenshot.

The Average Star Data dialog also has a button to generate a new PSF image computed from average parameters. Note that this has nothing to do with the average PSF image generated when you export a synthetic PSF, which is much more accurate and robust.

Exporting CSV Files

The comma separated value format (CSV) is a well-known, de facto standard to share data among applications. It is a rather simple format consisting of plain text data items separated with a delimiter, usually a comma character. The fact that virtually all spreadsheet and database management systems support CSV has encouraged us to implement it, as the best way to export fitted PSF parameters that can be further analyzed with other software applications.

CSV files exported by the DynamicPSF tool in PixInsight have the following properties and organization:

• Exported .csv files are plain ASCII text files.
• The data delimiter is a single comma character (ASCII code 44₁₀ = 2C_16).
• The decimal separator is a single period character (ASCII code 46₁₀ = 2E_16).
• The first line of an exported .csv file contains the names of the data items, namely: "ViewId,StarId,Channel,Function,B,A,cx,cy,sx,sy,FWHMx,FWHMy,unit,r,theta,beta,MAD"
• Each one of the rest of lines contains the data for a PSF fit. Note that the first item of each line is the identifier of the source view.

The block of text below shows an example CSV file generated with DynamicPSF:

ViewId,StarId,Channel,Function,B,A,cx,cy,sx,sy,FWHMx,FWHMy,unit,r,theta,beta,MAD
NGC1808_L,28,0,Moffat,0.010567,0.037276,995.48,1533.27,4.54,4.21,3.69,3.42,px,0.926,50.29,4.53,2.823e-03
NGC1808_L,7,0,Moffat,0.010599,0.027557,1427.22,579.95,4.88,4.31,3.76,3.32,px,0.884,58.90,5.01,2.823e-03
NGC1808_L,22,0,Moffat,0.010589,0.035163,1917.89,1228.73,5.05,4.53,3.79,3.39,px,0.896,52.51,5.27,3.120e-03
NGC1808_L,4,0,Moffat,0.010674,0.017232,1520.77,772.01,4.66,4.19,3.69,3.32,px,0.899,55.69,4.76,3.146e-03
NGC1808_L,30,0,Gaussian,0.011444,0.408243,672.41,706.42,1.70,1.70,4.01,4.01,px,1.000,0.00,0.00,3.206e-03
NGC1808_L,19,0,Moffat,0.010816,0.017442,1798.46,772.15,4.74,4.15,3.75,3.28,px,0.877,57.03,4.77,3.376e-03
NGC1808_L,31,0,Gaussian,0.011924,0.648031,1876.11,380.91,1.67,1.67,3.93,3.93,px,1.000,0.00,0.00,3.386e-03
NGC1808_L,25,0,Moffat,0.010521,0.012280,1220.03,611.32,4.43,4.05,3.65,3.34,px,0.914,55.97,4.42,3.567e-03
NGC1808_L,6,0,Moffat,0.010647,0.016519,1320.62,714.95,4.87,4.38,3.68,3.32,px,0.901,57.33,5.17,3.825e-03
NGC1808_L,32,0,Moffat,0.010511,0.013898,1297.07,424.24,4.66,4.66,3.62,3.62,px,1.000,0.00,4.94,4.928e-03
NGC1808_L,11,0,Moffat,0.010764,0.008830,1437.17,1169.46,4.67,3.98,3.79,3.23,px,0.851,55.29,4.55,5.594e-03