NGC 5189 from Gemini Observatory South:
Deconvolution and HDRWaveletTransform in PixInsight

By Juan Conejero (PTeam)
With original raw data from the Gemini Science Archive


Introduction

In this document we describe a complete processing work that focuses on two powerful tools in PixInsight: HDRWaveletTransform and Deconvolution. We show how both tools integrate in an advanced processing workflow applied to deep-sky astronomical images in PixInsight.

HDRWaveletTransform is our implementation of the High Dynamic Range Wavelet Transform algorithm devised by Vicent Peris, a member of the PixInsight Development Team. This tool is described in a specific tutorial that includes a step-by-step processing example. You may also be interested in looking at some examples of application to astronomical and diurnal images. HDRWaveletTransform is an extremely powerful tool to unveil hidden structures in images of objects possessing large ranges of brightness, as we'll demonstrate later in this document.

Deconvolution is our implementation of Richardson-Lucy and Van Cittert deconvolutions with wavelet-based regularization and deringing algorithms. Regularized deconvolution works by separating significant image structures from the noise at each deconvolution iteration. Significant structures are kept and the noise is discarded or attenuated. This allows for simultaneous deconvolution and noise reduction, which leads to robust deconvolution procedures that yield greatly improved results when compared to traditional or less sophisticated implementations. In-depth descriptions and references to the sources of the algorithms, with an example of application to high-resolution lunar images, can be found in a specific tutorial. The present work focuses on the Deconvolution tool applied to deep-sky images.

Both tools are state-of-the-art implementations of advanced image processing algorithms and techniques, representative of the power and innovation that are characteristic of the PixInsight platform.


The Image

For this example we have used an image of NGC 5189 acquired at the Gemini Observatory South with the GMOS instrument (Gemini Multi-Object Spectrograph). The image has been formed with original raw data downloaded from the Gemini Science Archive, and is used here with the kind permission from Gemini Observatory and NOAO.

On the Files and Links section you'll find data files and links to relevant resources related to this image, both on our website and on Gemini Observatory's website. For further information and comparison of obtained results, you may want to visit the corresponding entry for NGC 5189 on Gemini Observatory's image gallery.


NGC 5189 — Gemini South data processed in PixInsight with the tools and techniques described in this document.
Gemini Observatory / Association of Universities for Research in Astronomy.
Click on the image to see a full-resolution JPEG version.

Here is a brief description of this fascinating object, excerpted from Gemini's image release page: "NGC 5189, a chaotic-looking planetary nebula that lies about 550 parsecs (1,800 light-years) away in the southern hemisphere constellation Musca, is a parallelogram-shaped cloud of glowing gas. The GMOS image of this nebula shows long streamers of gas, glowing dust clouds, and cometary knots pointing away from the central star. Its unruly appearance suggests some extraordinary action at the heart of this planetary nebula."

We have used raw data taken with Hα, OIII and SII filters. To form a RGB composite image, we have mapped Sulphur to red, Hydrogen to red-orange, and Oxygen to cyan. We describe the color combination procedure in detail later on this document.

Special thanks to Vicent Peris (PTeam) for downloading and preparing the raw data for me, which has saved a lot of time.

Needless to say that processing this fantastic image has been a great pleasure. Working with such high-quality data is a privilege for anyone who loves astronomy and image processing.


Files and Links

In this section you can download several files and follow relevant links related to the image and data used in the tutorial. The PSM file and the individual raw FITS images can be used to reproduce the entire processing described on this document. All FITS files are 32-bit integer images compressed as .bz2 archives.


Acknowledgement

Based on observations obtained at the Gemini Observatory, which is operated by the Association of Universities for Research in Astronomy, Inc., under a cooperative agreement with the NSF on behalf of the Gemini partnership: the National Science Foundation (United States), the Science and Technology Facilities Council (United Kingdom), the National Research Council (Canada), CONICYT (Chile), the Australian Research Council (Australia), CNPq (Brazil) and SECYT (Argentina).


Step 1: Load the Raw Images

Our processing example begins loading the three raw images in PixInsight (Figure 1). Raw images are the Hα, OIII and SII components that we'll combine to obtain a starting RGB working image. Each raw image has been stored as a FITS file in 32-bit integer format. On Section 2 we have included a number of links that you can follow to download these files.

Each raw image is the result of averaging three raw data files that have been downloaded from the Gemini Science Archive. Individual frames have been extracted with the DS9 software package and have been aligned and integrated in PixInsight with the DynamicAlignment and PixelMath tools.


32-bit Integer Format Required

In this example, we are going to apply high dynamic range techniques, deconvolution algorithms and aggressive nonlinear transformations that involve very complex numerical procedures. With the 32-bit integer format we'll work in a very large numerical range of 232 discrete sample values, which are sufficient to perform such radical transformations. Another excellent choice would be the 64-bit IEEE floating point format, also supported by the PixInsight platform. The 32-bit floating point format provides less than 27 discrete values; it could be used, but a larger range is always preferable when it comes to perform very complex and long processing works.

Figure 1— The three raw FITS images loaded in PixInsight, corresponding to Hα, OIII and SII components, respectively.
(Click on the image to enlarge)


Defining Screen Transfer Functions

Since the raw images are linear, we can't see their contents directly on the screen. So the very first step is to open the ScreenTransferFunction interface, which will allow us to explore the images without changing them. A screen transfer function (STF) in PixInsight is a special histogram transform that is applied to pixel data just before they are rendered on the screen, but a STF does not modify actual image data in any way. STFs allow us to work with linear images in an easy and completely transparent way, just as we can do with nonlinear images.

On Figure 2 you can see the ScreenTransferFunction interface in action to define STF parameters for the Hα raw image (the rightmost image shown on the figure with the Ha identifier). Note the small midtones balance value (0.04), necessary to make all the image contents clearly visible. Since at this stage we are only interested in carrying out the initial alignment of the three image components, we haven't refined STF parameters any further. The STF shadows clipping value could be slightly increased to yield a much pleasant screen representation.

Figure 2— Fine tuning the screen transfer function for the Hα raw image with the ScreenTransferFunction interface.
(Click on the image to enlarge)

Now instead of repeating the same procedure for the remaining two images, we can tell PixInsight that the same STF parameters should be applied to all open image windows. This can be done by selecting the Screen Transfer Functions > Copy To All Views option from the image context menu, as shown on Figure 3.

Figure 3— Assigning the same STF to all images with the image context menu.


Step 2: Register the Raw Images

Now we'll open the DynamicAlignment (DA) interface to align the OIII image with respect to the Ha image (Process > ImageRegistration > DynamicAlignment, also available under the ImageRegistration category on the Process Explorer window). DA needs you to select a couple of images that will be mutually registered. You must click on a first image when the mouse cursor is a check mark with a "1" number. This first image will be the source image for DA, which will be used as the registration reference image. On Figure 4 you can see the appearance of the mouse cursor just before clicking on the Ha image of this example.

Figure 4— Selecting the source registration image with DynamicAlignment.

Then a second image must be selected (the mouse cursor changes to a check mark and "2"). This will be the target image for DA. At the end of the process, a duplicate of the target image will be translated, rotated, scaled or distorted, as necessary to register it over the source image. Note that the source and target images will not be modified by DA.

DynamicAlignment is a semi-manual image registration tool able to register images with subpixel accuracy and using one, two, or three and more alignment points. In the first two cases, an affine transform is used to register the target image over the source image. With just one alignment point, only translations are used (shifting the target image in both plane directions). With two alignment points, the target image can be translated, rotated, and scaled in a single direction. With three or more alignment points, DA uses a completely different numerical device, known as two-dimensional surface spline, or thin plate. A surface spline can adapt to arbitrary local distortions to perform an extremely accurate image registration work. The more alignment points defined, the more complex registration problems can be resolved, and with higher accuracy.

In all cases, DynamicAlignment is designed to use stars or other relatively small image features as alignment points. A centroid is calculated automatically for each alignment feature by fitting pixel sample values to a prescribed profiling function, which can be selected on the Sample Generation section of the DA interface. By default, a Gaussian function is used. In the present example, we'll use just default DA parameters, which usually work without flaws for most deep-sky images. DA can be used with linear raw images without problems, since it includes a special detection algorithm that automatically adapts to local brightness and contrast conditions.

In general, we recommend using at least three alignment points with DynamicAlignment. One or two points can be used in special circumstances, for example when there are no more valid points available, or to speed up iterative procedures applied to hundreds of images via scripts, when one knows in advance that no more precision is required. To maximize accuracy, the alignment stars must be selected to define a triangle covering a large portion of the image. Unsaturated, relatively small and round stars are the best candidates, although DA can work with elongated stars and other arbitrarily-shaped features.

On Figure 5 you can see the DynamicAlignment tool in action for the present example. We could use just three alignment points because there are no local distortions to correct for, since we are working with CCD images of the same object. However, we have preferred four alignment stars to enforce a precise correction for rotation. The chosen stars define a sort of rhomboid that extends well over both dimensions of the images. Significant differences up to 1.5 pixels have been found between the three raw images for some alignment stars.

Figure 5— Registering the SII raw image over the Hα image with DynamicAlignment. Note the four alignment stars.
(Click on the image to enlarge)

Note that we have created a process icon (align_raw_images on Figure 5) before clicking the Execute button on the DA interface. This has allowed us to reuse the same process instance to align the OIII image. The four alignment stars were automatically found and repositioned because the differences between the SII and OIII images are small.


Step 3: RGB Combination

The registered OIII and SII images have been assigned the OIII_registered and SII_registered identifiers, respectively, by DynamicAlignment. Before continuing, we'll close the original OIII and SII images (since they are no longer needed) and rename their registered counterparts to OIII and SII, respectively, for the sake of simplicity of the PixelMath expressions that we'll write. To change a view identifier, either double-click on its view selector (the vertical tab at the left edge of the image window), or select Image > Identifier.


Defining a Color Palette

For narrowband images, we have to define an arbitrary color palette to map individual filter components to RGB primaries. Indeed we are trying to represent wavelengths that are invisible to our eyes, and we often have two or more components at the same end of the spectrum, as Hα and SII in the present case, which complicates the situation.

The basic color palette that we have used for this example is as follows:

Color assignments must be carefully planned and implemented to optimize the representation of image structures, but achieving at the same time a plausible rendition from an aesthetical point of view. These are the color formulas that we have applied in this example:

To combine the three raw images as a RGB image, we'll implement the above formulas with the PixelMath tool.


RGB Combination with PixelMath

On Figure 6 you can see the PixelMath interface loaded with the corresponding expressions. Note that we have applied the PixelMath process globally to generate a new RGB image with the combination of raw components.

Figure 6— The PixelMath interface has been used to generate the combined RGB image from individual Hα, OIII and SII raw images.
(Click on the image to enlarge)

The PixelMath expressions are, respectively for the red, green and blue channels:

SII + 0.78*Ha
0.15*Ha + OIII
OIII

Note that we have added fractions of the SII and Ha raw images that sum more than one, so the result could be saturated, which would imply total loss of data. The same is true for the second expression (green channel). This happens, in fact, but just for the cores of a few bright stars that are already saturated in the individual components. This is because we are working with raw linear images, which are in general weak, and especially SII in this example. To verify that no critical regions are close to saturation, the following PixelMath expression can be used:

$target > 0.98

This expression is true (nonzero) for every pixel greater than 0.98 in the normalized [0,1] range (where 0 is black and 1 is white). For each color channel, nearly saturated areas will be pure white after applying this expression to a duplicate of the image, while the rest will be pure black.

Figure 7 shows the combined RGB image, after cropping and rotating it 90 degrees counter-clockwise, with a suitable screen transfer function applied. Note the strong red cast on the sky background; this is very easy to fix with the PixelMath process, as we'll see in subsequent steps.

Figure 7— Inspecting the combined RGB image with the help of the ScreenTransferFunction tool.
(Click on the image to enlarge)


Step 4: Linear RGB Working Space

In PixInsight, luminance/chrominance separations are performed into a colorimetrically-defined RGB working space (RGBWS) that can be specified independently for each image. These separations are of the utmost importance in many image processing tasks because the luminance supports visible image details, due to the characteristics and physiological limitations of the human vision system.

Before continuing, we must define a RGBWS to control how subsequent luminance/chrominance separations will be carried out. Two constraints must be taken into account for the image of this example:

Figure 8 shows the RGBWorkingSpace interface loaded with the parameters of the RGBWS that we have defined for this image. Note that we are free to define a triad of luminance coefficients that sum up more than one. The core PixInsight application automatically rescales them to form a unit vector when the instance is applied.

Figure 8— Defining a uniform linear RGB working space (RGBWS) with the RGBWorkingSpace tool.


Step 5: Background Neutralization

At this point we must inspect the histogram of the image before taking further decisions. You can see how it looks like with the HistogramTransform interface, on Figure 9. The image has a strong red cast in the shadows. This can be deduced from the relative displacements between the three histogram peaks. The red cast is of course clearly visible thanks to the STF that we are using so far —see Figures 7 and 10, for example.

A background neutralization procedure achieves a neutral gray color over sky background areas. For images whose histograms are dominated by pixels from free sky areas, this means that the three individual RGB channels will possess the same median values after neutralization. Graphically, this also means that the three histogram peaks will be aligned, since the median is a good estimate of the central peak's position in a uniform distribution. The image must be correctly calibrated and reasonably free from sky gradients due to light pollution and other additive effects that cause nonuniform illumination. Of course, in order to apply a direct background neutralization procedure, there must be some free sky regions represented on the image.

Background neutralization can be implemented in several ways. Some of them involve a manual adjustment work to graphically align the peaks of the RGB histograms. However, this is inaccurate and time-consuming. Here we'll implement an automatic procedure with the PixelMath tool in PixInsight.


Inspecting the Histograms

By inspecting the histogram of the image on Figure 9, you can see that the median of all blue sample values is the lowest of the three RGB channels, since its peak is the leftmost one. If necessary, accurate median values can be obtained with the Statistics tool, for example if the peaks are close together and not so prominent.

Figure 9— Inspecting the histograms with the HistogramTransform interface. The displacement of the red histogram peak toward the right denotes a strong red cast in the shadows (see also Figures 7 and 10). (Click on the image to enlarge)

Our goal is to achieve equal median values of background pixels for the red, green and blue channels. We'll force the red and green medians to be equal to the blue median, whose initial value is the lowest one. This will darken the background and, since our clipping points stay far away from the corresponding histogram peaks, no relevant data will be clipped but just some background noise.


Using Previews to Gather Sky Background Information

To ensure that we gather pure background pixel values, we'll define a preview over a free background region. For this example, we have created a small preview near the bottom right corner of the image; you can see it and its corresponding histograms on Figure 10. Note that on Figure 10 we have an active STF (refer to Figure 7), and that the histograms are much more peaked than those of the entire image (Figure 9). This is because we have a much more uniform distribution of pixel values on the preview, as expected for a region that is free from extended nonstellar objects and hence representative of the sky background in the image.

Figure 10— The background preview and its histograms.

Now the preview (Preview01 on Figure 10) must be extracted as an independent image. This can be done either by clicking on its view selector (the vertical tab at the left edge of the image window) and dragging it outside the view selector tray, or by selecting Preview > Make Image from the main menu. We'll change the identifier of the new image to background.


Implementation of Background Neutralization with PixelMath

The following PixelMath expressions will neutralize the background of the image:

$target - (Med( background ) - Med( background[2] ))
$target - (Med( background ) - Med( background[2] ))
$target

where background is, as we have said, the identifier of the image that we have created as a duplicate of Preview01. These expressions subtract the blue channel's mean sky level from both the red and green channels. They must be applied to the image with PixelMath's Rescale result option disabled. The result can be seen on Figure 11.

Figure 11— After background neutralization applied with PixelMath. (Click on the image to enlarge)

In cases where very high accuracy is required, or in difficult images that provide very few free sky background regions, one can define two or more previews and write PixelMath expressions that take the average of the medians from all of them. For example, with two previews:

$target - (Avg( Med( background1 ), Med( background2 ) ) - Avg( Med( background1[2] ), Med( background2[2] ))
$target - (Avg( Med( background1 ), Med( background2 ) ) - Avg( Med( background1[2] ), Med( background2[2] ))
$target

Note that the Avg() and Med() function calls are identified as invariant subexpressions by the PixelMath interpreter, so they are calculated only once prior to expression execution. At this point you may want to review a specific tutorial on PixelMath that includes a complete reference to PixelMath's syntax with examples.

If the initial histogram peaks happen to lie very close to the histogram's black point, then a sanity procedure is to add a small pedestal to neutralized pixels. This can be implemented in the PixelMath expressions to avoid any unwanted clipping; for example:

$target - (Med( background ) - Med( background[2] )) + 0.005
$target - (Med( background ) - Med( background[2] )) + 0.005
$target + 0.005

Figure 12— The combined RGB image after background normalization.
Since the image is still linear, a STF is necessary to see its contents without modifying it.

Now we have a combined RGB image with neutral background. Note that at this stage we still have a linear image, since the PixelMath expressions that we have applied in the previous step are linear transformations: we have just subtracted constant values from all pixels. Figure 12 suggests the true potential of this image, thanks to a suitable STF.


Step 6: Fixing Cosmetic Issues

The GMOS instrument is composed by three CCD sensors covering a 5.5 arcmin x 5.5 arcmin field of view. The junctures between adjacent sensors cause two small gaps. On our working image, these gaps are represented as several rows of dead pixels, resulting from the superposition of raw images. There are also many hot pixels throughout the entire image.

Before applying any edge enhancement procedure, as deconvolution in our case, it is advisable to fix these cosmetic problems. Otherwise they will cause severe ringing artifacts that are very hard to remove on the final processed image.

Figure 13— Dead pixel rows and hot pixels (cosmic rays?) in a crop of the combined RGB image, enlarged 4:1.


Fixing Dead Pixel Rows with a JavaScript Script

These rows of dead pixels (Figure 13) can be manually fixed with the CloneStamp tool. However, this is an extremely tedious task and it's difficult to do a good job without damaging the image, especially where the dead rows cross some nebular features. In PixInsight, the best way to fix these dead rows is by means of a simple script.

Below is the script that we have used for the present example. Basically, this script replaces bad rows with adjacent good rows, in a way quite similar to what an automatized clone stamp tool would do. Then it applies two iterations of a low-pass filter to smooth the fixed region.

/**
 * A script to fix bad pixel rows in the NGC 5189 Gemini image.
 * From the image processing tutorial "NGC 5189 from Gemini Observatory South:
 * Deconvolution and HDRWaveletTransform in PixInsight", by J. Conejero (PTeam).
 */

#include <pjsr/ImageOp.jsh>

function FixRows( image, rows )
{
   // For each channel of the image
   for ( var c = 0; c < 3; ++c )
   {
      // Select this channel.
      image.resetSelections();
      image.selectedChannel = c;

      // Central row in the set of bad rows.
      var y2 = (rows[c][0] + rows[c][1]) >> 1;

      // Select a good row above the first bad row, and copy it to the upper
      // half of the bad rows set.
      image.selectedRect = new Rect( 0, rows[c][0]-1, image.width, rows[c][0] );
      for ( var y = rows[c][0]; y <= y2; ++y )
      {
         image.selectedPoint = new Point( 0, y );
         image.apply( image );
      }

      // Select a good row below the last bad row, and copy it to the bottom
      // half of the bad rows set.
      image.selectedRect = new Rect( 0, rows[c][1]+1, image.width, rows[c][1]+2 );
      for ( var y = y2+1; y <= rows[c][1]; ++y )
      {
         image.selectedPoint = new Point( 0, y );
         image.apply( image );
      }

      // Convolve the bad rows region twice with a 5x5 Gaussian filter.
      image.selectedRect = new Rect( 0, rows[c][0]-1, image.width, rows[c][1]+2 );
      for ( var i = 0; i < 2; ++i )
         image.convolve( [1,  2,   3,  2,  1,
                          2,  7,  11,  7,  2,
                          3, 11,  17, 11,  3,
                          2,  7,  11,  7,  2,
                          1,  2,   3,  2,  1] );
   }
}

function main()
{
   // Get access to the currently active image window.
   var window = ImageWindow.activeWindow;

   // Do we have one?
   if ( window.isNull )
      throw Error( "There is no active image window!" );

   // Fix bad pixel rows by copying adjacent good pixel rows and convolving with a
   // low-pass filter. Each of the last three array arguments passed to the FixRows
   // routine specify the top and bottom rows of a rectangular region to be fixed.
   with ( window.mainView )
   {
      // Inform the core PixInsight application that we are about to modify the
      // image in this view. Without this, we have read-only access to the image.
      beginProcess();

      // Call the fix routine for each bad rows region.
      FixRows( image, [[ 65, 66], [ 65, 68], [ 65, 68]] );
      FixRows( image, [[747,747], [747,748], [747,748]] );
      FixRows( image, [[816,816], [816,817], [816,817]] );

      // Done modifying the image.
      endProcess();
   }
}

main();

To use this script, do the following:

The results of this script are not absolutely perfect, but they are in our opinion extremely good, as you can verify in Figure 14, and of course by trying it yourself with the image of this example. If you have the vice of perfection, you can fix a few remaining defects with the CloneStamp tool —we didn't do that, in fact, because the few remaining defects are of very minor importance.

Figure 14— Fixing dead pixel rows with JavaScript in PixInsight. A crop of the original image enlarged 4:1.

Mouseover:    [before applying the script]    [after applying the script]

The PixInsight platform provides a powerful environment for development of scripts in the JavaScript language. The PixInsight JavaScript Runtime (PJSR) includes a rich set of predefined objects and routines for image manipulation and processing. PJSR objects can be readily used and extended in script code to build extremely powerful and versatile tools. If you are interested in the scripting capabilities of PixInsight, the best learning resource is the official PixInsight JavaScript Runtime Documentation, and of course the development sections of PixInsight Forum.


Fixing Hot Pixels and with CloneStamp

In oversampled images, hot/cold pixels and cosmic rays are easy to identify because they are much smaller than the PSF. This means that every isolated group of, say one to four bright pixels, can be securely regarded as an artifact. This greatly facilitates fixing hot pixels with the CloneStamp tool, as happens with the image of our present example. In subsampled images, fixing hot pixels can be an extremely difficult to impossible task.

Figure 15— Fixing hot pixels with the CloneStamp tool.

The CloneStamp tool must always be used with care; it requires experience and good training. For this image we have used small brushes with radii between 1 and 5 pixels, and opacities ranging from 0.25 to 1.0. Low opacities are more accurate to fix artifacts over complex regions, for example over nebular regions. When hot pixels are located over delicate gradients, it is difficult to regenerate the gradients without causing problems that can be even worse than what we are trying to fix. When one feels that this is starting to happen, it is always better to leave things as they are.

To use CloneStamp, proceed as follows:

  1. Open the CloneStamp interface. It's under the Painting category in the Process Explorer window.
  2. Click on the image. This first click will inform CloneStamp that you want to start a new session for the specified target image.
  3. Ctrl+Click on a location over the image. This will define the origin, or source point, for all subsequent cloning actions.
  4. Click on a different location and drag to clone source pixels.
  5. Repeat steps 3 and 4 to continue cloning pixels with different source coordinates, if necessary.
  6. At any point, you can change brush parameters, which will apply to subsequent cloning actions. You can also undo, redo and delete actions with the local history control buttons.
  7. When you are done, please create a process icon on the workspace to save your work; this is very important because it will save you from having to repeat the whole process, in case it must be redone. Then click the green check mark button, or press Enter.
  8. If you forgot to create an icon at step 7, you still can create one from the History Explorer window: right-click on the image and select Load History Explorer, or press Ctrl+Alt+Shift+H; then select the CloneStamp history item in question and drag it to the workspace to create an icon.

In PixInsight, CloneStamp is a dynamic tool that records all cloning actions performed over an image. As happens with all processes in PixInsight, a CloneStamp instance can be saved as a process icon on the workspace. This allows us to replay all clone stamp actions at any point in our processing workflow. For example, suppose that after applying a complex and tedious CloneStamp session, you decide to undo several steps of the image's processing history. As long as you saved your work as a process icon, there is no need to repeat it: just drop the icon over the image, and all clone stamp actions will be replayed automatically for you. Or you can open the CloneStamp interface and drop the icon on its control bar: all clone stamp actions will be replayed as the initial part of a new CloneStamp session, which you can use to extend the clone stamp procedure with more actions.

Of course, a CloneStamp icon only works if the image hasn't been cropped since the icon was created or replaced. But even if the image has been cropped, it isn't difficult to write a script that modifies the recorded coordinates in a CloneStamp instance. As you can see, you always have full control over everything you do in PixInsight.


Step 7: Deconvolution

At this point we think that a brief description of deconvolution can be of interest to many readers. To understand how deconvolution works, we must start by describing how observed images originate from our instruments. Consider the image formation equation:

O = (I * P) + N

where:

The symbol '*' stands for the convolution operation. Convolution is a mathematical operation that causes a function to adopt the shape of another function. For example, if we convolve an image with a Gaussian function, all brightness variations in the image will be smeared according to the shape of a Gaussian distribution, which is very smooth. This has a blurring effect, hence the name Gaussian blur associated with that kind of convolution filters. Contrarily, if we convolve the image with a function that has a hard transition from negative to positive values, all brightness variations will be intensified locally, which causes an edge enhancement effect, also known as a sharpening filter.

So the image formation equation says that our observed image is just the result of a convolution with a point spread function (PSF). The PSF (P above) represents all instrumental imperfections and any fact or accident that might have affected the observing process, causing the observed image O to differ from the real image I. For ground-based astronomical observations, perturbations due to the Earth's atmosphere (seeing) form the greatest contributions to the PSF. Under ideal observing conditions —in absence of atmosphere, using a perfect telescope and a perfect tracking system—, the observed image would be diffraction-limited, which means that the PSF would be just the Airy disc corresponding to the instrument's aperture.

So, what is deconvolution? Ignoring the noise (N=0), we could write:

I = O / P

where, as you probably have guessed, the symbol '/' stands for the deconvolution operation. Deconvolution works by undoing the smearing effect caused to O by a previous convolution with P. In principle, if we'd happen to know the exact PSF, and there was no noise at all in the observed image, we could apply a direct deconvolution, as expressed above, to obtain something very close to the actual image. Of course, such a simple solution won't work in the real world, for three main reasons:

Deconvolution of real images is in general an ill-conditioned problem that requires sophisticated algorithms and implementations. Practical deconvolution algorithms include iterative and non-iterative solutions. Most common iterative deconvolution algorithms are Richardson-Lucy, Van Cittert and Maximum Entropy, and the most common non-iterative algorithm is Wiener deconvolution.

To obtain the PSF, there are basically two main strategies: estimate it by direct measurements on the image, for example by measuring the sizes and shapes of stars or other point-like sources, or use a blind deconvolution procedure. Blind deconvolution works by guessing different PSFs from actual image data and estimating the obtained improvement in the image after each try.

Finally, different regularization techniques can be used to prevent the adverse effects of the noise in the deconvolved solution, especially by controlling noise propagation and intensification in iterative deconvolution algorithms.



a

b

c


d

e

f


g

h

Figure 16— A deconvolution experiment in PixInsight. We simulate the image formation equation using a synthetic point source image, a Gaussian PSF, and uniform random noise. Then we deconvolve the simulated observed image (c) with several algorithms available in PixInsight (d, e, f, g, h).

a. Point source.
b. PSF.
c. Point source convolved with PSF + uniform noise.
d. Richardson-Lucy deconvolution, 50 iterations.
e. Regularized Richardson-Lucy, 100 iterations.
f. Wiener deconvolution, k-noise = 0.001.
g. Van Cittert, 6 iterations.
h. Regularized Van Cittert, 100 iterations.


Deconvolution Tools in PixInsight

As of writing this tutorial, the PixInsight platform offers two deconvolution tools, namely Deconvolution and WienerDeconvolution.

Deconvolution implements regularized Richardson-Lucy and Van Cittert deconvolution algorithms. Using wavelet-based multiscale analysis techniques, these algorithms separate significant image structures from the noise at each deconvolution iteration. Significant structures are kept and enhanced, while the noise is suppressed or attenuated. This allows for simultaneous deconvolution and noise reduction, which leads to robust, adaptive deconvolution procedures that in most cases outperform other implementations based on classical or less sophisticated design principles and techniques. In addition, our implementation includes a deringing algorithm that helps in fixing the effects of the Gibbs phenomenon, or ringing problems.

Our implementation of regularized deconvolution methods is based on many algorithms and ideas exposed by Jean-Luc Starck and Fionn Murtagh in Astronomical Image and Data Analysis (first and second editions, Springer, 2002 and 2006, resp.), with some additions necessary to write more versatile and adaptable tools, as well as slight modifications in the way noise is evaluated at each iteration.

The WienerDeconvolution tool has been written by our team member Carlos Milovic. Wiener deconvolution is a fast, one-step deconvolution algorithm in the Fourier domain. It is an interesting deconvolution tool that can be very useful to process any kind of images with relatively high signal-to-noise ratios. Wiener deconvolution is particularly suitable for lunar and planetary imaging. It is also good to shape star profiles, as correction for less-than-perfect stellar images due to minor optical aberrations or slight tracking errors, and to quickly find and test PSFs.

Both tools include selectable Gaussian, motion blur and custom PSFs, as well as a dynamic range extension feature that is very useful to fix saturated areas.

We strongly recommend reading a tutorial about deconvolution of high-resolution lunar images. That tutorial, which includes a complete, step-by-step processing example, covers many aspects of our algorithms and implementations in great detail.


Deconvolution Must Be Applied to Linear Images

Deconvolution actually only makes sense for linear images. Consider that no PSF (Point-Spread Function) can be valid simultaneously for all pixels of an image whose pixel values follow a nonlinear function of the intensities of the original sources. To better understand this, think of it in inverse terms: if we apply a nonlinear transformation to the raw image, the relative shapes and intensities of all image structures vary, so in theory we'd have to vary also the shape and size of the PSF as nonlinear functions of brightness for each pixel. This of course is not feasible.

There can be exceptions to this rule. Easy joke: but we can't remember one right now. Seriously, if you ever find yourself trying to deconvolve a nonlinear image, it might be that:


The Point Spread Function

As we have said, a good estimate of the PSF is crucial to perform a successful deconvolution. Basically, there are four methods to obtain a PSF:

  1. Direct measurements of star sizes on the acquired image.
  2. Accurate seeing measurements during image acquisition.
  3. Manual trial-error work.
  4. Theoretical models.

The preferred method is undoubtedly (1): when it comes to finding an accurate PSF, nothing can compete with direct measurements of what has been recorded on the image for point sources. When no stellar images are available, (2) can be used as a first approximation. The "manual" method (3) should be avoided, unless we have no other way to proceed; however, an experienced user may obtain a good approximation to the true PSF with patience and working systematically. Theoretical models can be useful in some cases; for example, the Airy disc corresponding to the observing instrument can be used as a first approximation to the PSF in lunar and planetary images acquired under good seeing conditions by integration of many short-exposure video frames.

The mean FWHM (Full Width at Half Maximum, see good descriptions here and here) is widely used as a standard measurement of image quality. For a given star in an image, the corresponding pixel values are fitted to a Gaussian function, from which the FWHM value is derived. The mean FWHM of the image is the mean of the same values calculated for a set of selected stars. The mean FWHM can be used to obtain an accurate PSF very easily, as we'll see now.

Fortunately, Gemini Science Archive data include mean FWHM values measured on the original raw images. For the NGC 5189 image of our example, we proceed as follows:

where σ is the standard deviation of the Gaussian distribution. Hence, assuming a Gaussian PSF, we have:

The standard deviation of 1.6 pixels characterizes the Gaussian distribution of the PSF that we'll use to deconvolve this NGC 5189 image.


The Choice of a Deconvolution Algorithm

For deep sky images, we generally use the regularized Richardson-Lucy (RRL) deconvolution algorithm (available in the Deconvolution tool). Of course, there's no reason for not exploring other options. As of writing this document, you can choose between regularized and non-regularized (classical) versions of the Richardson-Lucy and Van Cittert algorithms, along with Wiener deconvolution. From our experience, however, the best results are systematically obtained with RRL for deep sky images. For high-resolution lunar and planetary images, the regularized Van Cittert and Wiener deconvolution algorithms seem to be the best choices.


The Gibbs Phenomenon, aka Ringing

The Gibbs phenomenon poses a serious limitation on all edge enhancement image processing techniques. It is a hard problem, whose effects are particularly adverse for deep-sky images. It reveals as artifacts generated around high-contrast edges (jump discontinuities). For deconvolved stars and other small-scale bright structures, the Gibbs phenomenon causes ring artifacts that can easily lead to completely unusable results. Figure 17 shows an example of severe ringing. Note that the Gibbs phenomenon actually generates both dark and bright artifacts, respectively around bright and dark structures, although dark artifacts are generally the worst part of the problem. For strong discontinuities, as bright stars and other high-contrast features, ringing artifacts clearly reveal the oscillatory nature of the Gibbs phenomenon as alternating dark and bright rings.

Figure 17— A severe ringing problem caused after deconvolution. Note the strong ring caused by the hot pixels artifact, near the bottom edge of the crop. This clearly shows the importance of fixing all hot pixel artifacts before applying edge enhancement techniques.

Mouseover:    [Original crop]    [Classical Richardson-Lucy, 50 iterations]    [Wiener, k-noise=5×10-5]

This problem is particularly difficult to solve for linear images, where ringing can be a real nightmare. In the next sections, we'll show you our strategies to control ring artifacts caused by deconvolution of linear images. Of course, our solution is not perfect, but we think it is good enough as to allow practical deconvolution of most deep sky subjects. We are constantly thinking and working on this important problem, so you can expect to see significant improvements to our implementation in the medium term.


Building a Deringing Support

Our solution to the ringing problem is relatively straightforward: build a deringing support, which is basically a special star mask, and use it to command a deringing algorithm that works in tandem with deconvolution. A deringing support is similar to a mask that is nonzero only where ringing must be fixed. This comprises basically all of the stars in the image, along with other small-scale, relatively bright structures. For deep sky images, good deringing supports are vital to successful deconvolution. Despite being conceptually simple, building good deringing supports isn't usually an easy task because they are subjected to severe restrictions, necessary to prevent degradation of the deconvolution process. We'll see several examples in this section.

On Figure 18 you can see the support that we have used to deconvolve the NGC 5189 Gemini image of this example, along with the corresponding parameters of the StarMask tool that were applied to generate it. See the next section for a complete description of this important tool.

Figure 18— The deringing support for deconvolution, generated with the StarMask tool. (Click on the image to enlarge)

The complexity of a deringing support can vary widely, depending on the image being deconvolved. Basically, the main difficulties arise when there are small-scale, bright image structures that can cause ringing over nebular regions.


The StarMask Process

The StarMask process (MaskGeneration category) is the best tool in PixInsight to build star masks of many types, including deringing supports. StarMask operates by extracting the luminance from the target image (or a duplicate of the image if it's in the grayscale color space) and applying a multiscale algorithm that detects and extracts all image structures within a given range of scales (read sizes). The algorithm is based on the à trous wavelet transform and on a proprietary multiscale morphological transform.

It is very important that you understand the different parameters available to control generation of masks with the StarMask tool. Refer to Figure 18 to identify them on the StarMask interface.

Mode— You can select one between four available operation modes:


Star Mask


Structure Detection


Star Mask Overlay


Structure Detection Overlay

Figure 19— The four operation modes of the StarMask process.

Shadows, Midtones, Highlights— These parameters correspond to a histogram transform that is applied to the target image prior to structure detection and mask generation. In fact, this histogram transform is an important preparatory step in the StarMask algorithm. These parameters have default values of 0.0, 0.5 and 1.0, respectively, which define an identity transformation (no change). However, usually you'll need to apply lower values of the midtones balance parameter, especially working with linear images, mainly for two reasons:

On the example shown on Figure 20, note how the mask built with a midtones balance value of 0.1 includes much more stars and, at the same time, prevents inclusion of nonstellar structures from the brightest parts of the main object. We'll see later why this is particularly important for deconvolution deringing supports.

Increasing the Shadows parameter may also help to improve detection slightly; however, if you set it to an excessive value, clipping will occur in the shadows, which will prevent inclusion of dim structures. Generally, the highlights parameter is left with its default 1.0 value.

Figure 20— Two examples of the influence of the midtones balance parameter of StarMask.
In both cases the applied StarMask parameters are identical except midtones balance.

Target linear image (with active STF)

Star mask, midtones balance = 0.5 (default value).

Star mask, midtones balance = 0.1.

Scale— This parameter is the number of (dyadic) wavelet layers used to extract image structures. The larger value of Scale, the bigger structures will be included in the generated mask. Always try to set this parameter to the lowest value capable of extracting all required image structures; the range between 4 and 6 wavelet layers (scales from 16 to 64 pixels) covers virtually all deep-sky images.

Figure 21— The same mask generated with different values of the Scale parameter of StarMask.

Scale = 4

Scale = 5

Scale = 6

Growth— This is actually a category that comprises three StarMask parameters:

The first parameter (identified as Growth on the StarMask interface), controls the final sizes of all detected structures on the mask. The second and third parameters allow us to control the final sizes of small mask features with respect to larger features. See Figures 22 and 23 for some interesting examples.

Figure 22— The same mask generated with different values of the Growth parameter of StarMask.

Growth = 0

Growth = 1

Growth = 2

Growth = 3


Figure 23— The effects of different values of the Small-Scale Growth Compensation and Small-Scale Growth parameters of StarMask.

Scale = 5
Growth = 2
Small-Scale Growth Compensation = 0

Scale = 5
Growth = 2
Small-Scale Growth Compensation = 1
Small-Scale Growth = 4

Smoothness— This parameter determines the smoothness of all structures in the final mask. If generated with insufficient smoothness, the mask will probably cause edge artifacts due to abrupt transitions between protected and unprotected regions. Contrarily, excessive smoothness may degrade masking performance. In the case of a deringing support, finding a correct value for this parameter is very important. If in doubt, always prefer to exaggerate smoothness, because the effects of its lack are usually much worse.

Figure 24— The same mask generated with different values of the Smoothness parameter of StarMask.

Smoothness = 0

Smoothness = 10

Smoothness = 15

Smoothness = 20

Aggregate Structures— This parameter defines how individual image structures contribute to the mask construction process. This is a Boolean (enabled/disabled) parameter:

Figure 25— The effects of structure aggregation versus structure superposition in StarMask.

A mask generated by structure superposition.

The same mask generated by structure aggregation.

Binarize Structures—This parameter defines how the initial set of detected structures is truncated to differentiate the noise from significant structures. This is a Boolean (enabled/disabled) parameter:

Figure 26— Structure binarization versus structure normalization in StarMask.

A mask generated by structure normalization.

The same mask generated by structure binarization.

Threshold— This value is used to differentiate between noise and significant structures. Basically, all detected structures below this threshold will be considered noise and set to zero, and the rest will survive as significant structures. Obviously, higher thresholds will include less structures in the mask, and vice versa.

Figure 27— The same mask generated with different values of the Threshold parameter of StarMask.

Threshold = 0.2

Threshold = 0.1

Threshold = 0.05

Truncation— This value, in the range [0,1], is a highlights clipping point applied to the final mask (before multiplying by the Limit parameter, see below). It can be used to force the cores of bright structures to be pure white.

Figure 28— The Truncation parameter of StarMask.

Truncation = 1.0.

Truncation = 0.35.

Limit— This value, in the range [0,1], multiplies the whole mask after is has been completed, so it is useful to impose an upper limit for all mask pixels. Many deringing supports generated by structure binarization work better with lows limit values, between 0.1 and 0.5.

Figure 29— The Limit parameter of StarMask.

Limit = 1.0.

Limit = 0.25.


Deconvolving the NGC 5189 Gemini Image

Once we have generated a reasonably good deringing support with StarMask, it's time to test it with the Deconvolution tool. You know that your deringing support is "reasonably good" when you have gained experience with the practical usage of these tools, after having processed quite a bunch of images. We won't tell you lies: deconvolution isn't an easy task, and a good job usually requires a lot of testing, patience, and systematic work. Definitely, forget fast food solutions here, or you'll be wasting your time.

Figure 30 shows the Deconvolution tool in action for the NGC 5189 Gemini image of our example. Note the deringing support being used (the star_mask image) and the parameters on the Deconvolution interface. Note also that we make extensive use of relatively small previews to try out different deringing support images and deconvolution parameters. Previews are essential for this task, since deconvolution involves quite processor-intensive routines. A multicore or multiprocessor system is highly recommended; of course, the Deconvolution process will use all processors available. By default, Deconvolution creates high-priority threads, so don't expect your system to be very responsive during a long deconvolution procedure. If this behavior is a problem for you, it can be changed with the standard Preferences process (Edit > Global Preferences > Parallel Processing and Threads).

Figure 30— Deconvolving the NGC 5189 Gemini image with the Deconvolution tool in PixInsight. (Click on the image to enlarge)

On Figure 30, note that the standard deviation of the Gaussian PSF has been set to 1.6 pixels, which is the value that we derived from the Gemini GMOS plate scale and published FWHM values for this image. Note also that we have used the Regularized Richardson-Lucy deconvolution algorithm, which is our preference for deconvolution of deep-sky images. In what follows, we'll show you how to avoid some important pitfalls that usually arise during practical deconvolution. The rest is completely up to your experience, good taste, and common sense.


Apply the Correct Number of Deconvolution Iterations

Too few iterations won't improve your images too much. Too many iterations will lead to artifacts. When in doubt, it is always advisable to stop the process prematurely, since deconvolution artifacts are very objectionable. Applying an excessive number of iterations is a particularly dangerous risk with regularized deconvolution algorithms, since regularization is able to keep the noise perfectly controlled during hundreds of iterations.

The optimal number of iterations depends on the applied algorithm and on some characteristics of the image. As a general rule of thumb, a number between 30 and 50 seems to be a good choice for most images with regularized Richardson Lucy. Subsampled images usually deconvolve with just 10 or 20 iterations. Oversampled images, as is the case with this NGC 5189 example, usually require more iterations. You should try with several values and compare the results carefully.

Figure 31— A crop of the image enlarged 3:1, deconvolved with the regularized Richardson-Lucy algorithm and different amounts of deconvolution iterations. The temptation of applying too many iterations is particularly dangerous with regularized deconvolution algorithms, since noise intensification can remain perfectly controlled during hundreds of iterations.

In this example, 50 iterations are sufficient to show everything without overprocessing the image. Always be alert to stop deconvolution before it leads to "filamentary" artifacts inside bright features. These artifacts are conspicuous with 250 iterations in this example.

Mouseover:    [Original]    [20 iterations (too few)]    [50 iterations (correct)]    [100 iterations (too many)]    [250 iterations (horror)]


Optimize the Deringing Support

Always keep in mind that deringing actually works by degrading deconvolution performance. The deringing support must be tuned in order to keep such degradation effect restricted to stellar images and other bright features that cause objectionable ringing artifacts.

In practice, we always must tolerate ringing to some degree. The Gibbs phenomenon not only occurs around stars; it actually occurs around any discontinuity in the image, including all nebular structures that we love to see improved by edge enhancement techniques. Ironically, edge enhancement in digital images works thanks to the ringing problem that we are trying to fight here!

You can't remove the effects of the Gibbs phenomenon completely in the deconvolved image. That isn't the purpose of good deringing techniques, but just controlling (not suppressing) the amount of ringing originated around stars and other features where ringing cannot be tolerated.

When too much deringing is applied on regions where no deringing is necessary, deconvolution unevitably leads to artifacts. As in the case of applying too many iterations, improper deringing causes bright deconvolved features to adopt filamentary shapes. The reason is that deringing works by limiting decreased pixel values at each deconvolution iteration. When excessive deringing blocks darkening of pixels around deconvolved bright features, brightened pixels tend to group into islands inside the bright features. Figure 32 shows a good example of this effect, which must always be avoided.

Figure 32— Excessive deringing tends to generate isolated islands of brightened pixels inside bright nebular features, where deringing isn't necessary or must be kept to a minimum. The deringing support must be carefully built to minimize this effect. In this example, the deringing support #3 (see below) is excessive and is causing unnecessary degradation of deconvolution performance.

Mouseover:

[Original crop]    [Deconvolved, no deringing]

Deconvolved results:
[With deringing support #1 (good)]    [With deringing support #2 (good)]    [With deringing support #3 (bad)]

Deringing supports:
[Deringing support #1]    [Deringing support #2]    [Deringing support #3]

To optimize deringing supports, the Limit parameter of the StarMask process can be used to impose an upper limit of pixel values for the whole deringing support image. The deringing support #2 in the example of Figure 32 has been generated using structure binarization and a Limit value of 0.25. For the deringing support #1, which gives excellent results, we used structure normalization and a Limit value of 0.8.


Optimize Wavelet Regularization

As we have explained earlier in this document, regularized deconvolution applies a sophisticated, wavelet-based noise reduction algorithm at each deconvolution iteration. This allows us to deconvolve only significant structures in the image, avoiding intensification of noise, which is a major drawback of other implementations of deconvolution.

For the regularization mechanism to work properly, the user must define appropriate values for a set of regularization parameters. Basically, you must define:

Here are a few guidelines to properly handle these parameters for deep-sky images:

Figure 33— Wavelet regularization must be optimized to avoid noise intensification with no degradation of deconvolution for small image features. In this example, the third mouseover link below shows an excellent result, where virtually no noise increment can be observed after deconvolution. The fourth mouseover link shows an oversmoothed result due to excessive regularization.

Mouseover:

[Original crop]    [Deconvolved, no regularization]    [Deconvolved, good regularization]    [Deconvolved, excessive regularization]

[Good regularization parameters]    [Bad regularization parameters]


Refine Your Result Systematically

Don't expect to achieve a well deconvolved image by just a few clicks. Deconvolution of deep-sky images is a complex task that requires a lot of patience and methodic work. Some images can be relatively easy, but many others can be really challenging.

Basically, the amount of time required to achieve an optimal deconvolved result depends on the number of mistakes that you make. As you gain experience and training, you'll make less mistakes, your evaluations will be more accurate, and your results will be better and require less hours of hard work.

In PixInsight, it is very important that you use previews systematically. Create several previews over interesting parts of the image and use them to try out different sets of parameters. For example, a good parctice is to define a preview, process it, and make an image with it ( Preview > Make Image). Then you can place the image side to side with the preview and compare it with the result after changing some parameters. Another way to do this is to define a preview and clone it several times (Preview > Clone Preview); by cycling them you can perform a sort of bliking comparison that can be very useful.


Step 8: Chrominance Enhancement with ATrousWaveletTransform

Now that we have deconvolved the luminance of the image, how about the chrominance? On one hand, you may have noticed that some image features, particularly stars, have lost part of their chrominance. The reason is obvious: since deconvolution has increased contrast locally, the luminance is stronger for small-scale features where they have been brightened, and a stronger luminance implies less color saturation. On the other hand, at this point the chrominance of this image seems to need an overall boost, in our opinion.

There are many ways to improve the chrominance; the most obvious is through a color saturation curve using the CurvesTransform or ColorSaturation tools in PixInsight. We could do that with good results, but here we want to try something conceptually more elegant.

Figure 34— Applying a large-scale hyperbolic transfer function with the ATrousWaveletTransform tool. Note the overall improvement in the chrominance. (Click on the image to enlarge)

The ATrousWaveletTransform (ATW) process in PixInsight includes a large-scale transfer function feature that can be used in a similar way to the well-known DDP procedure. However, our implementation has an important advantage over DDP: it allows you to define the precise scale at which the hyperbolic function is to be applied. This means that, unlike DDP, ATW's large-scale transfer function does not affect small-scale image structures, but just those structures that are represented in the residual layer of the wavelet transform.

On Figure 34 you can see how we have applied a hyperbolic curve to the 16-pixel wavelet scale of the 5189 Gemini image, and on Figure 35 you can evaluate the outstanding results. Why does this work so well for the chrominance at this point? Mainly for two reasons:

Figure 35— Improving the chrominance of the deconvolved image with a large-scale hyperbolic transformation applied with the ATrousWaveletTransform process in PixInsight. Note the dramatic improvement in color saturation.

Mouseover:

[Before the large-scale hyperbolic function]    [After the large-scale hyperbolic function]


Step 9: Nonlinear Transformation

In order to exploit this image as the original data deserve, an aggressive nonlinear transformation is necessary to show the dimmest parts of the object. For example, look at the delicate nebular features that fade out smoothly, as feathers. All of these structures are very weak and definitely need a good push.

On Figure 36 you can see the applied HistogramTransform instance and the resulting image. Not surprisingly, the shadows clipping and midtones balance values are similar to those that we have been using so far as part of a screen transfer function (STF) that has allowed us to work with the linear image so far —see for example Figure 12. The small clipping applied at the shadows (0.0083) does not invalidate the background neutralization that we have carried out before (see Step 5), since it is being applied to the RGB/K combined channel of HistogramTransform, and hence equally to the three channels of the image.

Figure 36— Nonlinear transformation with the HistogramTransform interface working in real-time preview mode. (Click on the image to enlarge)


Step 10: High Dynamic Range Wavelet Transform

Now that we have a nonlinear image, where we can see all of those subtle structures that make an absolutely wonderful image, let's see what hides behind the fog. This sentence may seem intriguing at first sight, but in a few moments you'll fully understand what we mean with "fog".


The high dynamic range problem

Now that we have almost finished our processing work, it's time to face a complex problem that arises quite frequently with CCD images of bright deep-sky subjects: in the linear raw image, there are many different numerical values that support significant image structures over very bright and very dark regions. While we are not dealing with a real high dynamic range image, since we are not integrating multiple exposures of the same subject, we indeed have a problem related to the intrinsic dynamic range of the raw data. This is what we call a high dynamic range problem.

In practical terms, a high dynamic range problem means that we cannot reproduce all relevant image structures on the available display devices and media, such as computer monitors and printers. Most output devices work with eight bits per color channel at most. This means that we have to render the image using just 256 gray levels per color channel, which are clearly insufficient. In addition, output devices usually implement color spaces that cover just a small fraction of the theoretical RGB working spaces used to perform complex image processing tasks. To solve a dynamic range problem means to achieve a suitable compression of the existing range of numerical values in the image, so that we can render appropriately the whole image on low dynamic range devices.

There are different strategies and solutions to solve dynamic range problems. In the days of film and chemical darkrooms, astrophotographers devised clever masking methods for contrast reduction in prints of bright deep-sky objects (see for example B. D. Wallis & R. Provin, A Manual of Advanced Celestial Photography, Cambridge, 1988, pp. 189-197). Specific digital techniques, as the well-known DDP method (Digital Development Process) developed by Kunihiko Okano, have been widely implemented for dynamic range compression.


The HDRWaveletTransform Process

While simple techniques like DDP may help in relatively easy cases, we definitely prefer more sophisticated and refined algorithms and implementations. One of them is the High Dynamic Range Wavelet Transform algorithm (HDRWT) devised by Vicent Peris, a member of the PixInsight Development Team. The implementation of this algorithm in PixInsight is the HDRWaveletTransform process, which we'll use in the present example. At this point you may be interested in reading a description of HDRWaveletTransform with a practical example, and also in a gallery of example images processed with this innovative tool.

Let's explore this fantastic tool without more preambles. Figure 37 shows the applied parameters. Now you know what unveils when the fog disappears.

Figure 37— After HDRWaveletTransform.
(Click on the image to enlarge)

It is important to point out that the HDRWaveletTransform process must be tried out on previews that include the full range of brightness values that is present on the whole image; otherwise the obtained results on the preview won't be identical to what will be achieved for the image after applying the same instance. In practice, this means that unless we have a relatively small object surrounded by large free sky regions, the recommended procedure is to define previews occupying the entire image. This is a limitation of the HDRWT algorithm, and is a consequence of its multiscale nature, which has quite complex implications. In fact, as of writing this document, this is the only process that is not strictly previewable in PixInsight.


Number of Wavelet Layers

The most important and critical parameter of the HDRWaveletTransform process is the number of layers. This parameter refers to the number of contiguous wavelet layers to which the HDRWT algorithm is to be applied, and must be greater than one. The more layers, the larger structures will be preserved (read protected) by the algorithm. On Figure 38 you can see a comparison of results obtained from seven to four layers. The choice depends on what structures one wants to give more relevance, and of course some good taste and common sense are also important here. Our choice has been six wavelet layers (Fig. 38c). Five layers (Fig. 38d) is also a good option, but we think that with six layers the algorithm yields a more balanced rendition of the object. Less than five layers (Fig. 38e) give too flat results, in our opinion.

Figure 38— HDRWaveletTransform: testing several values of the number of layers parameter.


38a— Before HDRWaveletTransform.


38b— HDRWaveletTransform, 7 wavelet layers.


38c— HDRWaveletTransform, 6 wavelet layers (our choice).


38d— HDRWaveletTransform, 5 wavelet layers.


38e— HDRWaveletTransform, 4 wavelet layers.


Wavelet Scaling Function

The rest of HDRWaveletTransform parameters have been left with their default values, except the scaling function, which has been changed to 3x3 Gaussian. This is the wavelet scaling function used internally by HDRWT to perform wavelet transforms with the à trous algorithm. Peaked or sharp scaling functions are good at resolving small structures, while smooth or flat scaling functions work better at larger scales. Varying the scaling function is a good way to fine tune the effect of the HDRWT algorithm. A comparison can be seen on Figure 39.

Figure 39— HDRWaveletTransform: testing different wavelet scaling functions. In all cases the HDRWT transformation has been applied to six wavelet layers. The rest of parameters with their default values.

Mouseover:

[3x3 Gaussian (our choice)]    [3x3 Linear Interpolation]    [5x5 B3 Spline]

For your reference, this is the 3x3 Gaussian kernel used internally by HDRWaveletTransform:

0.04 0.20 0.04
0.20 1.00 0.20
0.04 0.20 0.04


Step 11: Fixing Doughnut Saturated Stars with PixelMath

A few bright stars have acquired a ring-shaped profile reminiscent of a doughnut. This has only happened inside bright stars that are saturated in the original linear raw data. You see an example on Figure 40, corresponding to a bright star near the center of the image.

These ring artifacts are a side effect of subtracting large-scale structures from saturated data. Fortunately, they can be fixed very easily by recovering the original brightness profiles just for these saturated stars, with the help of a special mask. We'll use the PixelMath process to carry out this task.

Figure 40— A doughnut saturated star. The star's profile has been transformed into an inverted, ring-shaped curve.

First, we need a mask that is white just for the doughnuts. Since we know that the doughnuts correspond to stars that are saturated in the raw images —but not in the processed image—, we can profit from the fact that the linear raw images are darker than the processed image everywhere except for these saturated features. With this in mind, here is the PixelMath expression that builds the mask:

raw[0] > $target[0] && raw[1] > $target[1] && raw[2] > $target[2]

where the raw identifier corresponds to the RGB combined image, and $target is meant to be the current processed image. This expression is true (nonzero) for every pixel of the RGB combined image that is brighter than its counterpart in the processed image, for all channels. Figure 41 shows the employed setup.

Figure 41— Building the doughnut mask with PixelMath.
(Click on the image to enlarge)

To make the doughnut mask more efficient, we blurred it slightly by removing its two first wavelet layers with the ATrousWaveletTransform process, and then activated it for the processed image before applying the following PixelMath expression:

raw

This expression replaces all processed pixels with (brighter) raw pixels only where the doughnut mask is not black. See the result on Figure 42.

Figure 42— Fixing the doughnut stars (ring-shaped saturated stars) with a masked PixelMath instance.

Mouseover:

[Before fixing doughnuts]    [Doughnuts fixed with masked PixelMath]


Step 12: Background Noise Reduction With ACDNR

The last step in our processing work is a simple noise reduction applied just to low-SNR regions of the image, as the sky background and the outer regions of the objects. This is very easy to implement with the ACDNR (Adaptive Contrast-Driven Noise Reduction) process in PixInsight.

First we must define a luminance mask for the ACDNR algorithm. The mask is very restrictive and blocks almost the whole object. It can be seen on Figure 43.

Figure 43— Building the ACDNR luminance mask.
(Click on the image to enlarge)

With this mask, we apply a simple ACDNR instance to the luminance and chrominance of the image. The applied parameters can be inspected with the corresponding icon of the PSM file that accompanies this tutorial. Figure 44 shows a comparison where you can evaluate the result of this final noise reduction procedure.

Figure 44— Background noise reduction with ACDNR (enlarged 4:1).

Mouseover:

[Before noise reduction]    [After noise reduction with ACDNR]