This post outlines the processing I used for my February 2013 gallery photograph. Processing consisted of the following steps, all implemented with MATLAB software I wrote except for the registration and interpolation step.
The spatial denoising step reduces detector spatial noise (i.e., fixed pattern noise) by dark subtraction and flat fielding using the formula
where s' is the output denoised subframe, s is the input noisy subframe, f is the flat field master and d is the dark subtraction master.
Registration and interpolation
I use PixInsight's StarAlignment process to register and interpolate the dithered, denoised subframes. This step is represented by the formula
where s' is the output registered and interpolated subframe, s is the input unregistered subframe and r is the registration reference subframe. I use the Bicubic Spline interpolation method with a clamping factor set to 0.3 to minimize undershoot artifacts (i.e., ringing) .
Rejection and integration
The rejection and integration step robustly combines registered and interpolated subframes for signal-to-noise ratio improvement and rejection of spurious noise structures. This step is represented by the formula
where i' is the output integration, S is the set of input registered subframes, r is the integration reference subframe, σl and σh are respectively the low and high sigma rejection thresholds and ρl and ρh are respectively the low and high range rejection thresholds.
Input pixels s are initially rejected by range thresholding as follows
Range thresholding removes zero pixels along the edges of the subframes introduced for alignment purposes by the registration and interpolation step.
After range rejection, input pixels s are normalized as follows
where s' is the normalized pixel value, δ(r) and δ(s) are respectively the Rousseeuw/Croux Sn scale estimator of the reference and input subframe and ξ(r) and ξ(s) are respectively the median location estimator of the reference and input subframe.
Normalized pixels s' are rejected by pixel stack sigma thresholding as follows
where ξ and δ are respectively the median and mean absolute deviation about the median of the pixel stack. Sigma thresholding is repeated on each pixel stack until convergence is achieved.
Sigma rejection removes spurious noise structures such as cosmic rays and airplane and satellite trails.
The weighted mean of the normalized, accepted pixels across all pixels stacks S' forms the output subframe as follows
where i' is the output integration and W are weights. The weight of pixels from a given subframe s' is as follows
where δ(s') and σ(s') are respectively the mean absolute deviation about the median and an estimate of the background Gaussian noise of the unregistered, uninterpolated counterpart of subframe s'. The counterpart rather than the subframe itself is used to avoid biased noise estimates due to interpolation.
The deblurring step attempts to compensate for the low-pass filtering performed by the telescope optics and the registration and interpolation step. This step is represented by the forumla
where i' is the output deblurred image, i is the input blurry image, α is an estimate of detector gain, σ is an estimate of standard deviation of detector Gaussian noise, β is a detail coefficient bias, and τ1 and τ2 are detail coefficient thresholds.
The deblurring step applies a function to the wavelet detail coefficients of an unormalized Haar wavelet decomposition of the input image. The function, empirically derived, classifies wavelet detail coefficients in three categories based on their magnitude relative to a mixed Gaussian-Poisson noise estimate: small, medium and large. Classification is defined by the τ1 and τ2 thresholds. Small and large coefficient values remain invariant. Medium coefficient values are scaled by the value 1 + β / 2^(j - 1) where 0 < β and 0 < j is the wavelet refinement level. Gaussian exponential functions are used to smooth the transition between categories as follows
where dj' is the output detail coefficient on level j, dj is the input detail coefficient on level j and cj is the input coarse coefficient on level j.
For the February 2013 photograph I set β, τ1 and τ2 equal to 0.7, 6 and 600, respectively. I set α equal to 10% of detector gain to roughly account for subframe averaging. I set σ equal to an estimate of the background Gaussian noise of the input image. See my blog post here for more information on background noise estimation.
The tone mapping step applies a midtones transfer function and a linear blend of several contrast limited adaptive histogram equalizations with differing radiuses to the input image as follows
The histogram equalization is masked with a luminance mask to reduce noise amplification in the shadow regions of the image. My implementation of histogram equalization employs several improvements to reduce histogram binning quantization artifacts.
 F. Luisier, "The SURE-LET Approach to Image Denoising", Swiss Federal Institute of Technology Lausanne, EPFL Thesis no. 4566 (210), 232 p., January 8, 2010.
 K. Zuiderveld, "Contrast Limited Adaptive Histogram Equalization", in P. Heckbert, Graphics Gems IV, Academic Press, 1994.
 P. Rousseeuw, C. Croux, "Alternatives to the Median Absolute Deviation", Journal of the American Statistical Association, 88(424):1273-1283, 1993 December.