The Region Around NGC 7000 and IC 5070:
ATrousWaveletTransfrom and HDRWaveletTransform in PixInsight

By Juan Conejero (PTeam)
With original raw data acquired by Thomas W. Earle


Introduction

This example focuses on the practical usage of the ATrousWaveletTransform (ATWT) and HDRWaveletTransform (HDRWT) tools in PixInsight. We show how these tools can be applied to accomplish the tasks of detail enhancement and high dynamic range compression, respectively, to process a wide field deep sky CCD image.

The HDRWT tool is our implementation of the High Dynamic Range Wavelet Transform algorithm due by Vicent Peris, a member of the PixInsight Development Team. HDRWT allows us to manage the dynamic range of astronomical images with unprecedented power and accuracy.

ATWT is a sophisticated process built around the à trous discrete wavelet transform algorithm (JL. Starck and F. Murtagh, Astronomical Image and Data Analysis, Springer, 1st ed., 2002, p. 29). This tool can be used to perform a wide variety of image restoration and noise reduction tasks.

We have described both tools and their practical usage in a number of processing examples and tutorials. In the Files and Links section you'll find links that you can follow to read those works, along with a number of interesting resources on wavelets and their applications to image processing.

In addition to these tools, we show some interesting techniques for image integration and evaluation of noise and signal-to-noise ratios, implemented as scripts in the JavaScript language.


The Image

Thomas W. Earle has kindly given us permission to use one of his CCD images for this processing example. It is a wide field image of the region around the North America and Pelican nebulas, NGC 7000 and IC 5070, respectively.


The Region Around the North America and Pelican Nebulas— raw data acquired by Thomas W. Earle and processed in PixInsight with the tools and techniques described in this document. Click on the image to download a 2500x2500 JPEG version (4.8 MB).

Image data:

Individual raw images were calibrated in Maxim DL. After calibration, they were registered, integrated and processed applying the tools and techniques described in this document. The entire processing work was carried out with the PixInsight platform running on a Linux workstation (Fedora Core 6 with the KDE desktop environment).

I want to express my gratitude to Thomas W. Earle for allowing me to work with his fantastic data.


Files and Links

In this section you can download several files and follow relevant links related to the image and the data used in this document. The PSM file can be used to inspect the exact parameters of every applied process. We have also included some important resources for those who want to learn more about wavelets and their applications to digital image processing.

Figure 1— The three raw FITS images loaded in PixInsight, with an active screen transfer function. (Click on the image to enlarge)


Step 1: Image Registration

We begin loading the three raw images in PixInsight (Figure 1). Calibrated raw images have been stored as FITS files in 16-bit integer format. With 4096×4096 pixels —the original sensor size—, each of them occupies 32 megabytes.

Before continuing, a suitable screen transfer function (STF) must be defined and applied to the three images. The characteristics and roles of STFs in PixInsight have already been described in preceding examples and tutorials. A STF is necessary for visual inspection of linear images on the screen, but it does not modify image data in any way.

The very first processing step that we must carry out is image registration. There are differences in position and rotation between the three raw images —including a slight drifting applied intentionally between individual exposures to minimize the contributions of small sensor defects— that must be corrected. Images must be translated and rotated as necessary to achieve their mutual superposition with subpixel accuracy. Once registered, we'll be able to integrate them into a single image with improved signal-to-noise ratio, as we'll see in the next step.

On Figure 2 you can see the DynamicAlignment tool being used to register image #2 over image #1. When three or more images are being registered for subsequent integration, it is very important to register them with respect to a unique reference image. In this example, we have registered images #2 and #3 over image #1. This avoids cumulative errors that would occur in a chained registration sequence; e.g. #2 over #1 and #3 over #2.

Figure 2— Image registration with the DynamicAlignment tool. Note the five registration stars covering the whole image, and the modified STF used to improve visibility of stars and graphics. The Ex and Ey values on the DynamicAlignment interface are the errors in pixels for the calculated centroid with respect to the predicted position. These values should be watched to evaluate the accuracy of the registration solution.
(Click on the image to enlarge)

For the images of this example, we have defined five registration points (or samples, as they are called on the DynamicAlignment interface). With three or more points, DynamicAlignment is able to correct for arbitrary translation, rotation, scaling and distortion. The more registration points, the more complex registration problems can be solved. For this example, three registration points would suffice: rotations between images are very small, and there is no scale change, since the three images have been acquired with the same instrument. For a wide field CCD image, our only concern regarding distortions might be differential atmospheric refraction, but this is not the case in this example.

On Figure 2, note that the five registration stars have been chosen so that they form a square covering basically the entire image. This is advisable, when possible, because in this way the registration algorithms use a more accurate framework. Note also that we have modified the active STF to darken the representations of the images. This greatly helps in identifying registration stars and improves visibility of vector graphics drawn by the DynamicAlignment tool.

In Figure 3 one of the registration stars is being inspected in more detail. The DynamicAlignment interface includes an automatic prediction system that can be used to evaluate the accuracy of the registration solution. The Ex and Ey values are the errors in pixels for the current registration point with respect to the predicted centroid position. These values are updated each time you define a new registration point, or when you move an existing one slightly with the mouse. The predicted position is calculated using all points defined except the point being predicted, and uses the same algorithms that are applied to perform the actual registration work.

Registration of the third image was a matter of a few clicks, as we had saved the DynamicAlignment instance as a process icon. Since the three images are very close in position (there are maximum differences of about three pixels), DynamicAlignment found the five registration stars automatically on the third image.

Figure 3— Verifying one of the registration stars and the registration error with respect to the predicted centroid position. (Click on the image to enlarge)


Step 2: Image Integration

Now we have three working images: the initial image #1, which has been our registration reference, and the registered images #2 and #3. The original (unregistered) images #2 and #3 are no longer necessary. For the sake of simplicity, we have changed the image identifiers to ngc7000_1, ngc7000_2 and ngc7000_3, respectively.

Image integration procedures allow to obtain a single image formed by combination of the data from two or more images of the same subject. The primary purpose of integration is to increase the signal-to-noise ratio in the combined image. Basically, the signal is constant throughout all of the images being integrated (with possible differences in intensity due to varying exposure times). However, the noise has a random distribution and hence tends to cancel out when two or more images are integrated.

In general, the best integration method, in terms of signal-to-noise ratio increment, consists of calculating the average of a set of images. The theoretical increment in the signal-to-noise ratio achieved by averaging N images is equal to the square root of N. In our case, it would be 1.732. On Figure 4 you can see the applied PixelMath instance and the resulting averaged image.

Figure 4— Average integration with PixelMath.
(Click on the image to enlarge)


Evaluation of the Increment in the Signal-To-Noise Ratio with JavaScript

We'll describe a procedure to evaluate the increment in the signal-to-noise ratio achieved in the integrated image, with respect to the original raw images. To this purpose, we'll apply a method for the calculation of the standard deviation of the noise in an image, assuming a Gaussian distribution of the noise.

Of course, the best way to implement these calculations in PixInsight is by means of a script written in the JavaScript language. Here is our implementation based on an iterative algorithm due to Jean-Luc Starck (see the reference in the script source code).

/**
 * Evaluation of Gaussian noise by the iterative k-sigma thresholding method.
 *
 * Reference:
 *    J.L. Stack, F. Murtagh, Astronomical Image and Data Analysis, Springer,
 *    1st ed., 2002, pp. 37-38.
 */

// Working array to gather noise at each iteration.
var N = new Array;

/**
 * Image callback routine to extract noise wavelet coefficients.
 *
 * s     is the value of the current sample, which is actually a wavelet
 *       coefficient in the first wavelet layer
 *
 * ks    k*sigma threshold, with k=3 and sigma is the standard deviation of the
 *       noise in the previous iteration.
 */
function GetNoise( x, y, c, s, ks )
{
   if ( Math.abs( s ) < ks )
      N.push( s );
   else
      this.currentSample = 1.0; // tag this coefficient as significant
}

/**
 * Iterative evaluation of the standard deviation of Gaussian noise.
 * We use k=3 and return an estimate to within 1% accuracy.
 */
function Noise( img )
{
   // Perform a one-layer wavelet transform using a linear scaling function.
   var W = img.aTrousWaveletTransform( [ 1/16, 1/8, 1/16,
                                         1/8,  1/4, 1/8,
                                         1/16, 1/8, 1/16 ], 1 );
   for ( var i = 0, s0; ; )
   {
      // Standard deviation of the noise for this iteration.
      var s;
      if ( ++i == 1 )
         s = s0 = W[0].stdDev(); // first iteration
      else
      {
         // Standard deviation of the set of thresholded coefficients.
         s = Math.stddev( N );

         // Return an estimate to within 1% accuracy.
         // As a safeguard, we allow no more than 20 iterations.
         var converges = (s0 - s)/s0 < 0.01;
         var enough = i == 20;
         if ( converges || enough )
         {
            console.writeln( (converges ? "Convergence reached" :
                                          "*** Warning: No convergence"),
                             " after ", i, " iterations." );
            return s/0.800; // back to the image space
         }

         s0 = s;
      }

      // Extract the set of noise pixels for the next iteration.
      N.length = 0;
      W[0].forEachSample( GetNoise, 3*s );
      if ( N.length < 2 )
         return 0;
   }
}

/**
 * Script entry point.
 */
function main()
{
   // Get access to the current active image window.
   var window = ImageWindow.activeWindow;
   if ( window.isNull )
      throw new Error( "No active image" );

   console.show();
   console.writeln( "<b>" + window.currentView.fullId + "</b>" );
   console.writeln( "Calculating noise standard deviation..." );
   console.flush();

   // Obtain an estimate of the noise standard deviation for the current view.
   var s = Noise( window.currentView.image );

   // Output results.
   console.writeln( format( "Noise standard deviation = %.8f, N = %u", s, N.length ) );
}

main();

This script estimates the Gaussian noise standard deviation for the current active view. To apply it, we have defined a preview over a region dominated by the sky background in the image. The same previews have been defined for the integrated image and for each original raw image (Figure 5). After execution of the script, these are the results:

Image Standard deviation of Gaussian noise SNR improvement in the integrated image
Integrated image 1.8228×10-4
Raw image #1 3.6755×10-4 2.02
Raw image #2 3.7079×10-4 2.03
Raw image #3 3.7287×10-4 2.05

The obtained mean increasing factor of 2.03 in the signal-to-noise ratio is well above the expected theoretical value of 1.73. The most probable cause for this difference is the fact that we have integrated one original image (#1) and two slightly translated and rotated images (#2 and #3), due to the necessary registration process. Pixel interpolation inevitably causes a slight loss of resolution, acting as a low-pass filtering procedure. A second cause could be the assumption of a Gaussian noise distribution, while the noise in CCD images normally follows a mixture of Gaussian and Poisson distributions, but the effect of this error would be negligible.

Figure 5— Previews defined for evaluation of Gaussian noise.
(Click on the image to enlarge)



Figure 6— To the left, a crop of the raw image #1. To the right, the same crop from the integrated image. Note the significant reduction in the visible amount of noise achieved by averaging the three raw images.


JavaScript Implementation of Image Integration with Pixel Rejection

Figure 4 shows a problem that unfortunately is quite frequent in images of moderate to long exposure times: airplane and satellite trails. In our example, we have three airplane trails crossing the image from side to side. Along with them, we have also several artifacts caused by hot CCD pixels and cosmic rays.

Removing all of these spurious data is relatively easy with a pixel rejection algorithm, as sigma-clipping for example, applied in tandem with the image integration procedure. We'll implement an image integration procedure that includes a simple but efficient pixel rejection algorithm, suitable to integrate small sets of images. Pixel rejection, however, exacts a price: the achieved signal-to-noise ratio is not as high as the ratio that can be obtained by averaging all pixels from the source images. Later we'll show you a technique to palliate this issue.

This is our implementation in JavaScript:

/**
 * A simple script to perform average image integration with asymmetric k-sigma
 * pixel rejection and generation of rejection map images.
 *
 * From the processing example:
 *    The Region Around NGC 7000 and IC 5070: ATrousWaveletTransfrom and
 *    HDRWaveletTransform in PixInsight, by J. Conejero.
 */

#include <pjsr/UndoFlag.jsh>

// Base identifier of the images to integrate
// The script will try to integrate images such as ngc7000_1, ngc7000_2, ...
#define BASE_ID      ngc7000

// Identifier of a preview to integrate. Leave it empty to integrate the
// whole images.
#define PREVIEW_ID   Preview01

// Number of source images to integrate.
#define IMAGE_COUNT  3

// Kappa value for rejection of bright pixels.
// ### You must fine tune this parameter for each specific integration/rejection problem.
#define KAPPA_BRT    0.12

// Kappa value for rejection of dark pixels.
// ### You must fine tune this parameter for each specific integration/rejection problem.
#define KAPPA_DRK    0.25

/**
 * Script entry point.
 */
function main()
{
   // Gather images and working dimensions.
   var width = 0;
   var height = 0;
   var images = new Array;
   for ( var i = 0; i < IMAGE_COUNT; ++i )
   {
      var id = #BASE_ID + '_' + (i + 1).toString();
      if ( (#PREVIEW_ID).length )
         id += "->" + #PREVIEW_ID;
      var view = View.viewById( id );
      if ( view.isNull )
         throw new Error( "No such image: " + id );

      images[i] = view.image;

      if ( i == 0 )
      {
         width = images[i].width;
         height = images[i].height;
      }
      else
      {
         // All images must have the same dimensions
         if ( images[i].width != width || images[i].height != height )
            throw new Error( "Incompatible image dimensions: " + id );
      }
   }

   // Integrated image in 32-bit integer sample format.
   var intWindow = new ImageWindow( width, height, 1, 32,
                              false, false, #BASE_ID+"_integration" );
   var intView = intWindow.mainView;
   intView.beginProcess( UndoFlag_NoSwapFile );
#define integration intView.image

   // Pixel rejection map image - bright pixels.
   var brtMapWindow = new ImageWindow( width, height, 1, 8,
                              false, false, #BASE_ID+"_rejection_map_bright" );
   var brtMapView = brtMapWindow.mainView;
   brtMapView.beginProcess( UndoFlag_NoSwapFile );
#define brtMap brtMapView.image

   // Pixel rejection map image - dark pixels.
   var drkMapWindow = new ImageWindow( width, height, 1, 8,
                              false, false, #BASE_ID+"_rejection_map_dark" );
   var drkMapView = drkMapWindow.mainView;
   drkMapView.beginProcess( UndoFlag_NoSwapFile );
#define drkMap drkMapView.image

   // Enable and initialize status monitoring.
   integration.statusEnabled = true;
   integration.initializeStatus( "Image integration", width*height );

   // Allow our users to abort the operation.
   console.abortEnabled = true;

   try
   {
      // For each row
      for ( var y = 0; y < height; ++y )
      {
         // For each column
         for ( var x = 0; x < width; ++x )
         {
            // Gather the stack of source pixels for the  current coordinates.
            var srcPixels = new Array;
            for ( var i = 0; i < IMAGE_COUNT; ++i )
               srcPixels[i] = images[i].sample( x, y );

            // Calculate the median of the source stack. This is an excellent
            // estimate of the central peak's position on the histogram.
            var m0 = Math.median( srcPixels );

            // Perform pixel rejection
            if ( m0 != 0 )
            {
               // Prepare to gather unclipped pixels.
               var pixels = new Array;
               var brtCount = 0;
               var drkCount = 0;

               // For each source pixel
               for ( var i = 0; i < IMAGE_COUNT; ++i )
               {
                  // Relative distance of this pixel to the central peak
                  var d = (srcPixels[i] - m0)/m0;

                  // Reject or accept this pixel
                  if ( d < 0 )
                  {
                     // This is a dark pixel
                     if ( d > -KAPPA_DRK )
                        pixels.push( srcPixels[i] );
                     else
                        ++drkCount;
                  }
                  else
                  {
                     // This is a bright pixel
                     if ( d < KAPPA_BRT )
                        pixels.push( srcPixels[i] );
                     else
                        ++brtCount;
                  }
               }

               // Update rejection maps.
               brtMap.setSample( brtCount/IMAGE_COUNT, x, y );
               drkMap.setSample( drkCount/IMAGE_COUNT, x, y );

               // Integrate the set of surviving (unclipped) pixels, or use
               // the median of all source pixels if all pixels were clipped.
               integration.setSample( pixels.length ? Math.avg( pixels ) : m0, x, y );
            }
            else
            {
               // In case we have a null (all zeros) pixel stack
               brtMap.setSample( 1, x, y );
               drkMap.setSample( 1, x, y );
               integration.setSample( 0, x, y );
            }
         }

         // Invoke the garbage collector after processing a whole row of pixels.
         // A good idea if we are integrating big images.
         gc();

         // Update the status monitor.
         integration.advanceStatus( width );
      }

      // Done with target views.
      drkMapView.endProcess();
      brtMapView.endProcess();
      intView.endProcess();

      // Show them.
      intWindow.show();
      brtMapWindow.show();
      drkMapWindow.show();
   }

   catch ( x )
   {
      // Kill them.
      drkMapView.cancelProcess();
      brtMapView.cancelProcess();
      intView.cancelProcess();
      delete intWindow;
      delete brtMapWindow;
      delete drkMapWindow;

      throw x;
   }
}

main();

This script is suitable to integrate relatively small sets of raw images; say no more than eight or ten. This is mostly because it implements a one-step pixel rejection algorithm. For larger sets, iterative rejection routines can yield superior results. Another limitation is that the script only integrates the first nominal channels of the source images, hence it expects to work with grayscale images. Interested readers with some programming skills can easily modify the script to overcome both limitations.

Along with the integrated image, the script generates two rejection maps for bright and dark pixels, respectively. The rejection maps are very useful to evaluate the quality of pixel rejection. Each pixel in a rejection map has a value proportional to the number of rejected pixels from all source images. If a rejection map pixel is zero, this means that all source pixels have been integrated (averaged) for the corresponding coordinates. Contrarily, If a rejection map pixel is one, it means that all source pixels have been rejected (this obviously can only happen if all source images are zero at the corresponding coordinates).

The most important parameters that must be fine tuned in the script are the KAPPA_BRT and KAPPA_DRK values. These are, respectively, the rejection thresholds for bright and dark pixels. A pixel is considered bright (dark) if its value is greater (less) than the median of all source pixels for its coordinates. The lower a threshold, the more restrictively performs the pixel rejection algorithm, and hence more pixels will be rejected (and the final signal-to-noise ratio will be inferior, too). A balance must be found so that spurious data are well rejected, as far as possible without rejecting too many random noise pixels. Figure 7 shows the obtained rejection maps and integrated image for our example. The optimum threshold values have been found to be 0.12 and 0.25, respectively for bright and dark pixels —remember that you must fine tune them carefully if you decide to try this script for your own images.

Figure 7— Results of the image integration script applied to a preview of the example NGC 7000 image.


7a— Pixel rejection map, bright pixels.


7b— Pixel rejection map, dark pixels.


7c— Integrated image.

In the resulting integrated image (Figure 7c), the weakest trail has not been completely removed, but it has been considerably mitigated. The rest of trails and artifacts have virtually disappeared. See Figure 8 for a detailed comparison.

Figure 8— Comparison between several integration procedures.

Mouseover:  [Raw image #3]  [our integration script]  [median integration]  [average integration]

In Figure 8, if you compare the results of a median integration and our integration script, you'll see little apparent differences. Besides the fact that our script has rejected airplane trails slightly better, the true difference lies in the signal-to-noise ratio. Using our Gaussian noise estimation script for high-SNR and low-SNR crops, eloquent results have been tabulated in Figures 9 and 10.

Figure 9— Noise evaluations for a crop of the image dominated by strong signal, after combination of the three raw images using several image integration methods.

Image Standard deviation
of Gaussian noise		 SNR Improvement
Our script 2.6804×10-4 1.89
Median integration 2.9972×10-4 1.69
Average integration 2.6601×10-4 1.90
Raw image #1 5.0589×10-4 1.00



Figure 10— The same test performed in Figure 9, but for a low-SNR crop.

Image Standard deviation
of Gaussian noise		 SNR Improvement
Our script 2.1001×10-4 1.72
Median integration 2.1589×10-4 1.68
Average integration 1.8165×10-4 1.99
Raw image #1 3.6217×10-4 1.00


Increasing the SNR for Rejected Pixels

Figure 9 shows that our integration script performs much better than a median combination for regions of the image dominated by strong signal levels (high-SNR regions). This is a logical result, since there are very few rejected pixels over those regions (have a look at the rejection maps shown in Figures 7a and 7b), so the script basically applies an average integration. However, Figure 10 says that this is not true for low-SNR regions, e.g. for the sky background, where our script performs basically like a median integration.

We'll show now a somewhat "tricky" technique that can fix this situation efficiently. If correctly implemented, this technique yields the same SNR improvement as an average integration, while keeping intact the rejection properties of our integration script.

The basic idea behind this technique is pretty straightforward: apply a low-pass filter exclusively to all rejected pixels in the integrated image. This will uniformize the values of all "black holes" that have been caused by pixel rejection. To implement this technique we must build a mask using both pixel rejection maps and the integrated image. We proceed as follows:

ngc7000_rejection_map_bright || ngc7000_rejection_map_dark

Note that the expression above applies a logical OR operation between both map images. The result of this operation is true (one) if any of both images is nonzero, false (zero) otherwise. The resulting mask can be seen on Figure 11.

Figure 11— Using PixelMath to combine both pixel rejection maps into a binary mask.
(Click on the image to enlarge)

$target * ~MTF( 0.003, Max( 0.004, blurred_duplicate ) )

The call to MTF (midtones transfer function) in the expression above is necessary because we are working with a linear image, which must be stretched. The 0.003 balance value was found after two or three tries. The call to Max avoids the contributions of abnormally dark pixels, like black rows and columns at the edges of the blurred_duplicate image. This expression must be applied to the result of the previous PixelMath instance (see Figure 11). Since we haven't changed its identifier, the mask image still is ngc7000_rejection_map_bright. The resulting mask is shown on Figure 12.

Figure 12— The mask after multiplication by the inverse of the stretched blurred_duplicate.
(Click on the image to enlarge)

With this procedure we have increased the signal-to-noise ratio for all pixels that were rejected by our integration script, since we basically have applied a low-pass filter to them. Removing the first wavelet layer from an image is equivalent to applying a low-pass filter, since the first wavelet layer supports small-scale —or high-frequency— image structures. The amazing result of this simple but efficient technique can be verified in Figure 13. For the same crop shown on Figure 10, we have improved the SNR by a factor of 2.08, clearly above the theoretical result of averaging integration for three images (equal to the square root of three, or 1.73). This indicates that our SNR improvement technique actually behaves like a selective noise reduction for low-SNR regions.

Figure 13— The same test performed in Figure 10, this time including the result of our special SNR improvement technique.

Image Standard deviation
of Gaussian noise		 SNR Improvement
Our script 2.1001×10-4 1.72
Our script with SNR improvement technique 1.7423×10-4 2.08
Median integration 2.1589×10-4 1.68
Average integration 1.8165×10-4 1.99
Raw image #1 3.6217×10-4 1.00

Mouseover:  [After SNR improvement]  [Before SNR improvement]


Step 3: Detail Enhancement with Wavelets

Deconvolution would be an obvious choice to continue processing the image of this example with PixInsight. The regularized Richardson-Lucy algorithm, available in our Deconvolution tool, can indeed perform very well for images like the present one, provided that a good estimate of the PSF and a suitable deringing support are used. For the best results, strongly undersampled images, as is the case in this example, can be upsampled prior to deconvolution. This involves working with upsampled PSFs. Among other reasons, this is important because a very small PSF cannot be accurately sampled when it is discretized to form a convolution kernel. 

But we pursue something completely different with this example. Instead of basing our processing on deconvolution, we'll work with the ATrousWaveletTransform tool, which will help us in achieving a strong enhancement of small-scale structures. Then we'll apply the HDRWaveletTransform tool to control the dynamic range of the image, in terms of overall and local contrast.

At this point, it is very convenient that you read a specific processing example on deconvolution and HDR wavelet transforms applied to deep sky images. Of particular relevance to the procedures that we are about to implement are the sections of that tutorial about the StarMask and HDRWaveletTransform tools.


Building a Star Mask

The Gibbs phenomenon and its important implications (ringing artifacts) have already been discussed in the document on deconvolution mentioned above; to understand why we need a star mask at this point, we encourage you to read the information given there. Summarizing, the Gibbs phenomenon manifests as artifacts generated around all discontinuities present in the image when we apply a high-pass filtering process. Without a robust methodology to avoid ringing problems, edge-enhancement techniques are useless due to the objectionable artifacts that they generate. 

Ringing is particularly problematic around stars and other small-scale, relatively bright image features. To overcome this problem, our "standard solution" consists of building and using special star masks, either to protect those small and bright image features, or to command special deringing algorithms (as the deringing feature of the Deconvolution tool). To build star masks, StarMask is a powerful and versatile tool in PixInsight. At the risk of being amazingly repetitive, we encourage you to read the mentioned processing example on deconvolution.

Figure 14— Building a star mask with the StarMask tool.
(Click on the image to enlarge)

Our star mask, and the applied StarMask parameters, are shown on Figure 14. We must point out that the StarMask process, due to particularities of the algorithms implemented, should be applied to main views (whole images), or to previews comprising the entire range of brightness present in the image. Otherwise, the results obtained on a partial preview may differ from the results obtained on the whole image with the same StarMask parameters. The safest way to try out StarMask is to define and use a preview covering the entire image. We haven't found, as of writing this document, a satisfactory workaround for this limitation.


Combined Luminance / Star Mask

The star mask will play a key role to avoid ringing artifacts, and also to prevent saturation of bright stars. However, this is not sufficient for the present example. Our intention is to boost small-scale features with an aggressive wavelet transform, as we'll see in the next section. This requires protecting low-SNR areas, as the sky background and all relatively dark regions, as well as all transitions between bright nebular features and the background. So what we need is a sort of multipurpose mask, which should be suitable to protect all stars and low-SNR regions at the same time.

Remember that a mask protects more where it is darker. So our mask should be darker were we need more protection:

The general procedure to build a combined mask that fulfills these requirements is as follows:

Figure 15— Building a luminance mask with the HistogramTransform tool working in real-time preview mode.
(Click on the image to enlarge)

mask * ~star_mask

On Figure 16 you can see how we have combined both masks in our example. Note on this figure that not all of the stars have been removed from the combined mask. Two facts are important here and must be taken into account:

Figure 16— Combining the luminance mask and the star mask with PixelMath.
(Click on the image to enlarge)


The Wavelet Transform

Now that we count with the protection of a good mask, the real fun begins with the ATrousWaveletTransform tool. Prepare to see a dramatic change in Figure 17.

Figure 17— Trying out ATrousWaveletTransform on a preview.
(Click on the image to enlarge)

The following facts must be pointed out for the ATrousWaveletTransform shown on Figure 17:

0.0625 0.125 0.0625
0.125 4 0.125
0.0625 0.125 0.0625
As you see, it's small and has a prominent central peak.

Noise reduction works on a per-layer basis in ATrousWaveletTransform, so you can define independent sets of noise reduction parameters for each wavelet layer. The kernel size parameter must be kept as low as possible to achieve the required noise reduction without damaging fine detail. In general, applying more filtering iterations and a small value of the amount parameter helps in preserving fine structures, and is preferable over a few iterations with a higher amount.

When applying noise reduction and layer biasing simultaneously for one or more wavelet layers, consider that you are working with two opposite forces, each one trying to increase what the other tries to decrease. It's as if you were applying a low-pass filter and a high-pass filter at the same time; you must tune each filter carefully to narrow the exact image structures that you want to enhance, keeping the noise controlled as much as possible.

Figure 18— The wavelet transform. What happens if we suppress some elements of the transform?

Mouseover:

[Processed with ATrousWaveletTransform with deringing + noise reductio + mask]

[Before the wavelet transform]     [No deringing]     [No noise reduction]     [Unmasked]

Figure 19— Wavelet noise reduction is better evaluated on a crop enlarged 2:1.

Mouseover:

[Processed with ATrousWaveletTransform with deringing + noise reductio + mask]

[Before the wavelet transform]     [No noise reduction]


Step 4: Nonlinear Transformation

So far the image is linear. Now it's time to apply a nonlinear transformation to make it representable on display media. The HistogramTransform tool must be used to perform this task in PixInsight. As you see on Figure 20, the applied transformation is quite aggressive (a midtones balance value of 0.006). These are the two histogram rules that you always must respect:

Figure 20— Nonlinear transformation with HistogramTransform.
(Click on the image to enlarge)


Step 5: High Dynamic Range Wavelet Transform

If conventional deep-sky images are just what you want to achieve, then you can bypass this step completely and go directly to noise reduction. Not too many people will tell you that your images are unnatural if you do so. But if you want to explore new paths to improve the expressive strength of your astrophotography, then stay with us and learn a few secrets about the dynamics of astronomical images. We prevent you that a lot of people will tell you that your images are overprocessed and unnatural. The choice is yours.

The HDRWaveletTransform tool implements a clever and extremely versatile algorithm due to our friend Vicent Peris. He knows much more than me about the dynamics and expressivity of images, so I'm acting here as a mere translator. Nevertheless, Figure 21 shows a few tries with this exciting tool. To know more about it, we think (this is the last time you read this here) that you should review another processing example on deconvolution and high dynamic range wavelet transforms, along with a specific example about HDRWaveletTransform.

Figure 21— HDRWaveletTransform applied at different scales. All process parameters have default values except number of layers (which has been varied as indicated below, and scaling function, which has been changed to "3x3 Linear Interpolation".

Mouseover:

[Before HDRWaveletTransform]

[HDRWT applied to 10 wavelet layers]     [9 layers]     [8 layers]     [7 layers]     [6 layers]     [5 layers]

My choice has been HDRWT applied to seven wavelet layers, but I also like the result of applying it to eight layers a lot, and I recognize that with eight layers we have a more "complete" image, in the sense that all structures are better represented. This is just a matter of personal preferences. With perhaps the exception of the 5 layers option (absence of large-scale features), the rest of possibilities seem perfectly valid in my opinion. Of course, nothing prevents you from combining two or more versions. Varying the scaling function is also a good way to fine tune the effects of HDRWaveletTransform, as illustrated on Figure 22.

Figure 22— Varying the scaling function of HDRWaveletTransform. In all cases the algorithm has been applied to seven wavelet layers.

Mouseover:

[HDRWT applied with 5x5 B3 Spline]     [3x3 Linear Interpolation]     [3x3 Gaussian]


Step 6: Noise Reduction

The last step of our processing example consists of a noise reduction applied with the ACDNR tool (Adaptive Contrast-Driven Noise Reduction). You have the employed setup on Figure 23, and a good comparisons on Figure 25. Note that ACDNR is being used with the luminance mask option enabled, which is essential to protect high-SNR areas and transitional regions that don't need as much noise reduction as the background. The mask is shown on Figure 24.

Figure 23— Noise reduction with ACDNR.
(Click on the image to enlarge)


Figure 24— Fine tuning the ACDNR mask with the real-time preview interface.
(Click on the image to enlarge)


Figure 25— Noise reduction with ACDNR.

Mouseover:

[After ACDNR]     [Before ACDNR]

Finally, after noise reduction the histogram is no longer dominated by the noise, which allows us to perform a final stretch to improve overall contrast. On Figure 26, note the three small peaks that are being clipped to the left of the main histogram peak. These peaks correspond to a few rows and columns of dead pixels that were generated at the edges of the raw images as a side effect of image registration.

Figure 26— Final stretch after noise reduction.
(Click on the image to enlarge)