Bayer Drizzle Issues?

[sharkmelley]

Are there any known issues with Bayer Drizzle?  I'm using it and finding no real difference in sharpness between the usual integrated image and the Bayer Drizzled image, which surprised me.  I'm using a Sony A7S on a Tak Epsilon (focal length 500mm) so there is plenty of potential for the Bayer Drizzle algorithm to make a difference.

I've used the both the BatchPreProcessing script and the BayerDrizzlePrep script and they give identical results but not perceptibly better than the usual integration.

So now I'm playing with some synthetic data, generating synthetic raws from a star filled image with added features such as single red, green and blue pixels, blocks of 2x2 red, green and blue etc.  In theory, just 4 of these synthetic raws (with single pixel offsets) should be sufficient for drizzle to provide an almost perfect replica of the original.  However I'm finding that in each of the 4 "subs", Drizzle is smearing each single coloured pixel into a block of 4.  This is not what I expected to see and might explain why Drizzled images look just as blurry as standard integrated images.  In fact the tests I've done so far indicate that the standard integration is actually doing a better job of reproducing my synthetic image than Drizzle. 

I can tidy up my 4 synthetic raws and make them available as a test example if that would be useful.

Mark

[sharkmelley]

So here's my test.  I've superimposed a pattern of coloured squares and stripes on a star background.  From this I created 4 synthetic raw files with single pixel offsets to cover a complete CFA 2x2 square.  These synthetic raws therefore sample every image pixel in every RGB channel.  The background stars serve to align the images during integration.  A standard integration and a bayer drizzle integration is then performed (scale=1 dropshrink=1.0).  If bayer drizzle is working properly it ought to be able to exactly reproduce the original image from which the synthetic raws were created.  In fact I've used both the bayer drizzle approaches provided by PI - one using the BatchPreprocessing script and one using the BayerDrizzlePrep script but both give very similar results so I present the BayerDrizzlePrep result below, magnified by 4 so the individual pixels are visible:



The standard integration (I used bilinear interpolation in the debayer) actually does a very creditable job.  This gives me confidence that my synthetic raws are correct.  The single coloured blocks and 2x2 coloured blocks come out well but slightly blurred.  The white stripes are resolved well and the green stripes look good.  The red and blue stripes fail to resolve because of the sparseness of the R&B pixels.

The bayer drizzle integration is actually worse than the standard integration.  The single coloured blocks now appear as 2x2 blocks.  The 2x2 coloured blocks no longer look like 2x2.  All the thin stripes now fail to resolve.

So am I doing something wrong?  Or is the PI bayer drizzle algorithm doing something different to what I would expect?

The test data can be found here:
https://drive.google.com/drive/folders/0B3Ky5pyZvsINeW42RjlHbkVLSDA

The files are all in PI XISF format to prevent the usual issues of "up-bottoming" or not with FITS files. The bayer pattern is "RGGB".  Frame1 - Frame4 are the "raw" files and I've provided synthetic master bias, dark, flat.  For convenience I've also provided the original from which the synthetic raws were produced.

Mark

[sharkmelley]

To clarify further.  The original image and the 4 synthetic raw files generated from it are 512x512 and all are full of stars.  Here is a stretched version:



The star background is what allows PI to register the images.  The added coloured dots and stripes are simply there as a device to give some detailed insight into what the bayer drizzle is doing to the data.

Mark

[Juan Conejero]

Hi Mark,

Thank you for your test. I am working on a new test with your data for a comparison between Bayer drizzle and standard image integration. I'll post the projects and scripts involved later, along with an analysis of results.

[Philippe B.]

Very strange
I made some nice images with Bayer drizzle and here are the comparison :
We can really see improvement in bayer drizzle version

[Juan Conejero]

Hi Mark,

I have completed an interesting test with your original image. Besides demonstrating that my implementation of drizzle is correct (and for that matter, of Bayer drizzle), the procedures and results of this test are very useful to understand how the drizzle algorithm works in practice, so I'm going to detail them here.

The test begins with a little script, whose source code can be seen here:

Code: [Select]

#define ORIGINAL_IMAGE_ID  "original"
#define NUMBER_OF_IMAGES   100
#define MAX_DELTA_PX       25
#define OUTPUT_BASE_DIR    "/home/juan/tmp/bayer-drizzle-test"
#define BAYER_PATTERN      "RGGB"

#include <pjsr/Interpolation.jsh>
#include <pjsr/UndoFlag.jsh>

function translateImageByIntegerPixels( image, dx, dy )
{
   image.shiftBy( Math.round( dx ), Math.round( dy ) );
}

function translateImageByFractionalPixels( image, dx, dy )
{
   image.translate( -dx, -dy );
}

function convertRGBImageToBayerRGB( view, bayerPattern )
{
   let Rx, Ry, Gx0, Gx1, Bx, By;
   switch ( bayerPattern )
   {
   // R G
   // G B
   default:
   case "RGGB":
      Rx  = 0; Ry  = 0;
      Gx0 = 1; Gx1 = 0;
      Bx  = 1; By  = 1;
      break;
   // B G
   // G R
   case "BGGR":
      Rx  = 1; Ry  = 1;
      Gx0 = 1; Gx1 = 0;
      Bx  = 0; By  = 0;
      break;
   // G R
   // B G
   case "GRBG":
      Rx  = 1; Ry  = 0;
      Gx0 = 0; Gx1 = 1;
      Bx  = 0; By  = 1;
      break;
   // G B
   // R G
   case "GBRG":
      Rx  = 0; Ry  = 1;
      Gx0 = 0; Gx1 = 1;
      Bx  = 1; By  = 0;
      break;
   }
   let P = new PixelMath;
   P.expression0 = "iif( x()%2 == " + Rx.toString() + " && y()%2 == " + Ry.toString() + ", $T, 0 )";
   P.expression1 = "iif( iif( y()%2 == 0, x()%2 == " + Gx0.toString() + ", x()%2 == " + Gx1.toString() + " ), $T, 0 )";
   P.expression2 = "iif( x()%2 == " + Bx.toString() + " && y()%2 == " + By.toString() + ", $T, 0 )";
   P.useSingleExpression = false;
   P.createNewImage = false;
   P.showNewImage = false;
   P.newImageWidth = 0;
   P.newImageHeight = 0;
   P.newImageColorSpace = PixelMath.prototype.RGB;
   P.newImageSampleFormat = PixelMath.prototype.f32;
   P.executeOn( view, false/*swapFile*/ );
}

function main()
{
   let originalWindow = ImageWindow.windowById( ORIGINAL_IMAGE_ID );

   for ( let i = 0; i < NUMBER_OF_IMAGES; ++i )
   {
      let window = new ImageWindow( 1, 1, 1, 16/*bitsPerSample*/, false/*floatSample*/ );
      let view = window.mainView;

      let dx = 2*MAX_DELTA_PX*Math.random();
      if ( dx > MAX_DELTA_PX )
         dx = MAX_DELTA_PX - dx;

      let dy = 2*MAX_DELTA_PX*Math.random();
      if ( dy > MAX_DELTA_PX )
         dy = MAX_DELTA_PX - dy;

      view.beginProcess( UndoFlag_NoSwapFile );
      view.image.assign( originalWindow.mainView.image );
      translateImageByIntegerPixels( view.image, dx, dy );
      view.endProcess();

      convertRGBImageToBayerRGB( view, BAYER_PATTERN );
      window.saveAs( OUTPUT_BASE_DIR + format( "/integer/bayer/test%04d.xisf", i+1 ),
                     false/*queryOptions*/, false/*allowMessages*/, true/*strict*/, false/*verifyOverwrite*/ );

      view.beginProcess( UndoFlag_NoSwapFile );
      view.image.assign( originalWindow.mainView.image );
      view.image.interpolation = Interpolation_Bilinear;
      translateImageByFractionalPixels( view.image, dx, dy );
      view.endProcess();

      convertRGBImageToBayerRGB( view, BAYER_PATTERN );
      window.saveAs( OUTPUT_BASE_DIR + format( "/fractional/bayer/test%04d.xisf", i+1 ),
                     false/*queryOptions*/, false/*allowMessages*/, true/*strict*/, false/*verifyOverwrite*/ );

      window.forceClose();

      console.noteln( format( "* %4d : dx=%+6.2lf dy=%+6.2lf", i+1, dx, dy ) );
   }
}

main();
This script generates two synthetic data sets from your original image. Both sets consist of copies of your original RGB image, translated randomly in both plane directions, and sampled as Bayer RGB data. The images in one of the sets are translated by integer distances in pixels, while the images in the other set are translated by fractional pixels. To run this script, first you must create two directories on your filesystem, namely:

OUTPUT_BASE_DIR/integer/bayer
OUTPUT_BASE_DIR/fractional/bayer

where OUTPUT_BASE_DIR must be defined appropriately in the script (see line 5). Run PixInsight, load your original image, and set its identifier equal to "original" (or alternatively, change the value of ORIGINAL_IMAGE_ID in line 2 of the script). Load the script in Script Editor and execute it. By default the script generates 100 images in each set. You can redefine this number in line 3, by changing the value of NUMBER_OF_IMAGES. Finally, you can also change the macros MAX_DELTA_PX, which is the maximum displacement applied in pixels (in absolute value), and BAYER_PATTERN to one of "RGGB", "BGGR", "GRBG", and "GBRG".

As I have said the script generates Bayer RGB images. With a simple change we could generate Bayer CFA monochrome images. Besides some disk space savings, the result would be exactly the same. I prefer Bayer RGB images because they are much easier to inspect visually.

After running the script, I have also generated fake master bias, dark and flat frames using the NewImage tool (bias and dark as pure black images of 512x512 pixels, and the master flat with a constant value of 0.5 and also 512x512 pixels). Then I have executed the BatchPreprocessing script on each data set. In BPP, I have selected the appropriate light frames (either on integer/bayer of fractional/bayer), the fake master calibration frames, and have enabled the CFA images and Bayer drizzle options. I have disabled image integration because I have performed this task manually after BPP. Here is a screenshot of BPP's dialog window just before running it:

After running BPP, we have a set of calibrated/debayered/registered images, a set of Bayer drizzle files, and a set of Bayer RGB frames, which are redundant in this case since we already have them after running our synthetic data generation script.

With these images, I have generated four test results:

1. ImageIntegration: Standard integrated image using synthetic frames with integer translations, debayered, and registered:

2. DrizzleIntegration: Bayer drizzle integrated image using synthetic frames with integer translations and drizzle files.

3. ImageIntegration: Standard integrated image using synthetic frames with fractional translations, debayered, and registered.

4. DrizzleIntegration: Bayer drizzle integrated image using synthetic frames with fractional translations and drizzle files.

Here is your original image for comparison:

In all cases the images are shown enlarged 6:1. Only your original image is shown with an active STF, since enabling it for the rest would make evaluation of results impossible. The result images are available for download here:

http://forum-images.pixinsight.com/20161019/drizzle/original.xisf
http://forum-images.pixinsight.com/20161019/drizzle/integration-integer.xisf
http://forum-images.pixinsight.com/20161019/drizzle/drizzle-integer.xisf
http://forum-images.pixinsight.com/20161019/drizzle/integration-fractional.xisf
http://forum-images.pixinsight.com/20161019/drizzle/drizzle-fractional.xisf

Here are screenshots of the ImageIntegration and DrizzleIntegration tools, as used during the test:

   

As you can see, for simplicity I have disabled image weighting and rejection in all test cases. This is a screenshot of the workspace with the original and test result images for a better comparison:

(click to see a full-size image)

For DrizzleIntegration I have used a 1:1 drizzle scale and no drop shrink. It would be interesting to repeat the same test with different scale and drop shrink parameters.

How does Bayer drizzle perform as a function of the number of input images? Here are some additional results with only the first 10, 20 and 50 synthetic frames, respectively:

Evaluation of Results

In (1) we see how a standard image integration has produced a reasonably good result for images translated by integer pixels. However, two problems become evident: color fringing (the white features are now greenish with strong color variations), and aliasing (the red and blue features at the top right corner of your artificial test pattern have not been resolved at the one-pixel scale).

In (3), which is the standard integration of images translated by fractional displacements, the result is so bad that it is, in fact, an excellent demonstration of why Bayer drizzle should always be used, when possible, instead of the traditional deBayer/registration/integration procedure. We can see strong aliasing and color fringing artifacts that simply destroy the original test pattern.

In (2), the result of Bayer drizzle for images translated by integer pixel displacements, there are no color fringing artifacts, as expected. However, and also as expected, the drizzle algorithm cannot recover structures at the one-pixel scale when the source images are not dithered by subpixel displacements. This is true because of the way drizzle works, and is true irrespective of the number of source frames.

Finally, the resulting image in (4), which is the Bayer drizzled image for frames translated by fractional pixel displacements, demonstrates that my implementation of the drizzle algorithm works well. It has worked remarkably well, in fact, to reconstruct your artificial test pattern with no color fringing artifacts and all structures quite correctly represented at the one-pixel scale.

This test is particularly good to show three important facts about drizzle:

- Drizzle requires correctly dithered images. Dithering must be performed by non-integer pixel displacements in both directions, ideally following a random pattern.

- Drizzle wants many images. The typical "the more, the better" recommendation is particularly true for drizzle.

- When the data set allows implementing it, the Bayer drizzle technique is always superior to a standard deBayer/register/integrate procedure for mosaiced data. The benefits of drizzle in this case are minimal color artifacts, absence of aliasing artifacts, and better reconstruction of small-scale image structures.

I want to point out an important fact about this test, which should not pass unnoticed for a complete evaluation of the results. In the data generation script I have used bilinear interpolation (see script lines 98 and 99) to translate the images by fractional pixel amounts. Actually, this does not reproduce the quality of a true data set acquired with real dithering, since pixel interpolation necessarily acts as a low-pass filtering process, degrading the data at small scales. I have used bilinear interpolation because it produces no ringing artifacts, but it is also one of the worst choices in terms of aliasing. With real high-quality data, the result in (4) should be significantly better without this limitation. We could try to overcome this problem in the test by not interpolating the data, which would require an inverse drizzling of the original image to generate the synthetic frames, but this can be somewhat complex to implement.

[vicent_peris]

Hi,

A very short definition of the drizzle algorithm would be "an image reconstruction technique based on fractional pixel displacements". Any discussion about drizzle should meet this definition.


Best regards,
Vicent.

[sharkmelley]

Hi Juan,

Thank you so much for the time and effort you have spent in your very thoughtful and very educational reply.  The script you wrote is very efficient for creating the synthetic raw frames - it took me a lot of effort without that script!

You certainly prove that for random fractional displacements the PixInsight Bayer Drizzle algorithm works well with a scale of 1.

However I disagree with your assertion that "the drizzle algorithm cannot recover structures at the one-pixel scale when the source images are not dithered by subpixel displacements."  I see no reason why bayer drizzle cannot work with integer displacements.  Indeed, I can prove this with my original 4 synthetic raws - I'm sure you will agree that without loss of generality, 4 carefully chosen synthetic raws are sufficient when using integer displacements. 

See the image below:



Deep Sky Stacker is able to take those 4 synthetic raws and reproduce all the structures in my test.  Furthermore, PixInsight also reproduces those structures if I set the drizzle scale to 2 instead of 1.

So my question is, for integer displacements, why does PixInsight not reproduce those structures at a scale of 1 when DSS is able to do so?  I now think I can see what is happening at a scale of 1.  What should happen is that we have a target image and each coloured pixel from the source images "rains down" onto the target image landing precisely on top of target image pixels.  However what PixInsight currently seems to be doing is that each source pixel hits the middle of a group of 2x2 pixels in the target image and thus becomes a block of 4.  Possibly there might be an unintentional extra offset of half a pixel for some reason. Hence the target image produced by PI is equivalent to the original image convolved with a 2x2 block.  If you look very carefully at the drizzled image produced by PI you will see this is true.  You can also drizzle integrate a single image of the 4 - you will see that everything is constructed from 2x2 blocks.



My guess is that a minor change to the implementation of the PI bayer drizzle algorithm will fix this issue for integer displacements and then PI would reproduce all the structures in the test image even with a scale of 1.  What I can't work out is whether or not that fix would also improve the quality of drizzled images having random fractional displacements.

In the meantime I think I will use PI Bayer Drizzle with a scale of 2.

Kind regards,

Mark

[Juan Conejero]

I see no reason why bayer drizzle cannot work with integer displacements.

It is not the Bayer drizzle technique, but the drizzle algorithm. In the first place, bear in mind that drizzle only makes sense for undersampled data. If all of the source images are mutually displaced by integer amounts without any rotation, source drops tend to fall over single output pixels entirely. If this happens, the algorithm cannot reconstruct image structures smaller than one pixel (undersampled structures), simply because it cannot gather enough information from the data set. In the case of Bayer mosaiced images, the Bayer pattern causes a lack of spatial resolution (besides possible undersampling caused by the optical system) due to sparse sampling, especially in the red and blue channels. Demosaicing algorithms try to recover the lacking pixels by interpolation. Bayer drizzle reconstructs the lacking pixels without interpolation from a sufficiently large data set with random fractional dithering, which allows it to gather data "between" source pixels optimally. You need to understand the drizzle algorithm well (and ideally attempt to implement it) in order to fully understand this essential requirement.

I'm sure you will agree that without loss of generality, 4 carefully chosen synthetic raws are sufficient when using integer displacements.

Not without loss of generality. Your 4 synthetic raws are indeed carefully chosen because they cover any Bayer pattern completely. You are working with an ad-hoc data set that is not representative of a general set of images mutually displaced by integer pixel amounts without rotation. You need random displacements to approximate a generalized data set. However, this is not an important fact in this discussion.

Possibly there might be an unintentional extra offset of half a pixel for some reason.

It is intentional. Fortunately, the ImageIntegration module, and hence the DrizzleIntegration tool, has been released as an open source product. See DrizzleIntegrationInstance.cpp, lines 347 to 350. Also lines 311 and 322 (although the latter two are actually irrelevant for this purpose). These offsets of half a pixel have to be added to the initial drop coordinates in order to be coherent with the coordinate system used by the StarAlignment process. SA represents local pixel coordinates with the origin at the center of the top left pixel of each image. This is necessary to support specular transformations (horizontal and vertical mirroring), where the images have to be mirrored with no pixel shift at the origin. Again, this is unimportant to this discussion, since it is just a convention about the origin of a coordinate system.

My guess is that a minor change to the implementation of the PI bayer drizzle algorithm will fix this issue for integer displacements and then PI would reproduce all the structures in the test image even with a scale of 1. What I can't work out is whether or not that fix would also improve the quality of drizzled images having random fractional displacements.

Your guess is correct, and I can work it for you: it destroys the implementation completely for fractional dithering. As a test I have recompiled the ImageIntegration module with lines lines 347 to 350 of DrizzleIntegrationInstance.cpp commented out, as well as comments to remove the 0.5 pixel offsets in lines 311 and 322 (although these two changes make no practical difference in this case, since the involved coordinates are only used with rejection enabled). Here is the drizzle integration result with the same 100 randomly, fractionally displaced images of the original test:

For comparison, this is the original result with the correct 0.5 pixel offsets in place:

And here is the result of the modified DrizzleIntegration tool with the 100 frames with integer displacements:

which reproduces the test pattern "perfectly", as you guessed.

These are the result images:

http://pixinsight.org/images/forum/20161019/drizzle/drizzle-fractional-wrong.xisf
http://pixinsight.org/images/forum/20161019/drizzle/drizzle-fractional.xisf
http://pixinsight.org/images/forum/20161019/drizzle/drizzle-integer-wrong.xisf

In the "drizzle-fractional-wrong" image, besides an evident degradation of the test pattern, the reconstructed real sky image is also worse. The "drizzle-fractional-wrong" image has an average PSF larger than the correct one by about a 4%, as measured with the FWHMEccentricity script. The reconstructed sky image with integer displacements, "drizzle-integer-wrong", is identical to the original image. This is the expected result, since as I said at the beginning of this post, with integer displacements each drizzle drop covers a single output image pixel, given that we are working with a scale of 1 and no drop shrink.

Certainly it would be possible to modify several processes, including DrizzleIntegration, to improve our drizzle implementation for integer displacements without rotation. However, that would be a wrong decision for several obvious reasons, including: drizzle works optimally with random fractional displacements, real-world images are almost never dithered by integer amounts (this may happen with high-precision mounts, but for dithering it has to be interpreted as a problem to solve, not as a desirable feature), and real-world images are normally subject to at least small rotations. This makes the current implementation, which works optimally for arbitrary displacements, rotations and local distortions, the most logical and efficient option.

In the meantime I think I will use PI Bayer Drizzle with a scale of 2.

That makes no sense in general, except if you have very large data sets of (optically) undersampled data and want to perform a drizzle x2. There is no reason to use Bayer drizzle with a scaling factor of 2 on a regular basis. If you prefer DeepSkyStacker's implementation, or if you think that DeepSkyStacker's implementation is better, then why not use DeepSkyStacker?

[sharkmelley]

Hi Juan,

Thanks once again for taking the time to write such a comprehensive reply and for running that experiment with the change in source code.  I wasn't aware that the code for that module was available - I'll certainly take a look at it to understand it better.

I totally agree that integer displacements are a special (perhaps pathological) case and that under normal acquisition conditions the displacements will be fractional.  So I agree with you that it is definitely more important that the code works optimally in the fractional case, especially in the presence of minor frame to frame minor rotations.  But I do find it unintuitive that this should preclude optimal behaviour in the case of integer unrotated displacements. 

I'll have to give that a lot more thought.

Mark

[vicent_peris]

I want to repeat two sentences that I consider of extreme importance:

In the case of Bayer mosaiced images, the Bayer pattern causes a lack of spatial resolution (besides possible undersampling caused by the optical system) due to sparse sampling, especially in the red and blue channels. Demosaicing algorithms try to recover the lacking pixels by interpolation. Bayer drizzle reconstructs the lacking pixels without interpolation from a sufficiently large data set with random fractional dithering, which allows it to gather data "between" source pixels optimally.

Drizzle works optimally with random fractional displacements, real-world images are almost never dithered by integer amounts (this may happen with high-precision mounts, but for dithering it has to be interpreted as a problem to solve, not as a desirable feature), and real-world images are normally subject to at least small rotations. This makes the current implementation, which works optimally for arbitrary displacements, rotations and local distortions, the most logical and efficient option.

Please note that there is no need to apply a bayer drizzle x2 since we are replacing the debayering by drizzle. And, please, never (and I say "never" very few times) use drizzle after debayering since you'll be applying drizzle to an already interpolated image (this does not make any sense!).

I think Juan is way too optimistic. Even with extremely good mounts you cannot have images moved by exact entire pixels (maybe if you shoot at very short focal lengths). This means that the test does not represent at all a real situation under the sky.


Best regards,
Vicent.

[sharkmelley]

Hi Juan,

I've taken a look at the code: https://github.com/PixInsight/PCL/blob/master/src/modules/processes/ImageIntegration/DrizzleIntegrationInstance.cpp

If I understand it correctly, it works as follows:
1) referenceRect defines a target pixel in the destination image
2) sourceP0-sourceP3 are its corners in source image coordinates
3) The code then identifies the pixels from the source image that will contribute into the destination pixel (lines 362-372).  These are the pixels that will be drizzled into the target pixel
4) The function GetAreaOfIntersectionOfQuadAndRect calculates the intersection of the drizzled pixel and the target pixel in order to weight the data.

For the examples below I assume that Drizzle Scale =1.0 and Drop Shrink = 1.0, for simplicity.

Consider the case where the target pixel is more or less coincident with the a source pixel (i.e. just like our integer offset example).  It should receive a full contribution (weight =1.0) from the source pixel.  However, the intersected area is calculated from sourceP0-sourceP3 which has been offset by 0.5.  As a result, the target pixel receives equal contributions (i.e. weighted by 0.25) from a group of 4 source pixels.  Equivalently, each source pixel contributes equally to a 2x2 block of target pixels.  This is the effect I noticed and commented on earlier.

But it is not just "integer offset" pixels that are adversely affected by the unintuitive weights.  Consider the case where the target pixel sits centrally on the intersection of a 2x2 block of source pixels. We would expect each of those 4 source pixels to contribute equally to the target pixel.  However, because the intersected area is calculated from sourceP0-sourceP3 which has been offset by 0.5, it instead receives a full contribution from just one source pixel.

Generalising the above argument and working through a few more examples it quickly becomes apparent that for every fractional offset, the weights are incorrect because they've been calculated from an area of intersection that has been offset by 0.5 pixels from the target pixel.

So at first sight, it seems to me that removing the offset of 0.5 is actually the correct thing to do.  However, you did try that and we saw an obvious degradation in resolved features in the random fractional offset test case.  That puzzles me and I don't have an answer for it yet.

It might be worth trying it out on real data (instead of the synthetic test) to see if removing the offset is generally a good or bad thing.  However I completely understand if you have already spent enough time experimenting with Bayer Drizzle 

Meanwhile, when I get time, I'll try to compile and run the code myself. 

Mark

[Juan Conejero]

Hi Mark,

After making many tests with different data sets and working conditions, and having analyzed the problem more in depth, I think you are right: the 0.5 offsets are incorrect. Thank you so much for your endurance and patience with me

Indeed, StarAlignment sets the origin of its local coordinate system at the geometric center of the top left pixel of each image, for the reasons I explained: this simplifies all geometric transformations, including specular transformations, since the pixel at x=0 y=0 does not rotate when the image is mirrored horizontally or vertically, or rotated around the origin. However, the coordinate system used by the rest of the platform, and specifically by all PCL image support classes, fixes the image coordinate system at the top left corner of the top left pixel. You can verify this graphically; just open an image, enlarge it well (zoom in to 20:1 for example), and see how the reported coordinates are xxx.5 and yyy.5 when you move the mouse to the center of any pixel.

My mistake was trying to compensate for this translation of 0.5 pixels for calculation of drop intersections. There is no need to compensate for this offset in drizzle, since the coordinates computed by the coordinate interpolation device (either a homography or surface splines) are valid irrespective of the location of the origin of coordinates, as the coordinate system does not change. So yes, the 0.5 offsets have to be removed. The relevant sections of the code should be changed as follows:

line 311:
               referencePixel.y = TruncInt( referenceRect.y0 /*+ 0.5*/ );

line 322:
                  referencePixel.x = TruncInt( referenceRect.x0 /*+ 0.5*/ );

lines 347-350:
//                sourceP0.MoveBy( 0.5 );
//                sourceP1.MoveBy( 0.5 );
//                sourceP2.MoveBy( 0.5 );
//                sourceP3.MoveBy( 0.5 );

line 390:
                           double value = m_data.source( TruncInt( p->x + 0.5 ), TruncInt( p->y + 0.5 ), c );

The last change, in line 390, is necessary because the two coordinate systems---from image registration and the platform default---are displaced by 0.5 pixels, as I have said above. So the source pixel to be addressed has to be corrected by the same offset, since here we are mixing two different coordinate systems.

This will be fixed in the next version 1.8.5 of PixInsight. Since it's going to take long (multiple issues with Qt 5.7), I'll release an update with this bug fixed next week. BTW, I plan on adding many improvements to my drizzle implementation, including direct support for Bayer CFA frames (no need to split into RGB CFA) and, much more importantly, new drizzle drop kernels. Right now it uses a fixed square kernel, but it will support arbitrary kernels including a user-defined PSF, Gaussian, Moffat and variable kurtosis functions.

Interestingly, this bug seems to have no practical repercussion for data with randomly distributed fractional dithering. That's why we haven't discovered it before. It has an obvious repercussion for data with integer displacements, as you have discovered.

So at first sight, it seems to me that removing the offset of 0.5 is actually the correct thing to do.  However, you did try that and we saw an obvious degradation in resolved features in the random fractional offset test case.  That puzzles me and I don't have an answer for it yet.

I am puzzled, too. In all of my subsequent tests with synthetic data and your test image, removing the offset does not degrade the performance of DrizzleIntegration for fractional dithering. The results are practically identical with and without the fix. See for example this test with the bug fixed and a newly generated synthetic data set:

In this test I have included also a Bayer drizzle x2 result. With 100 perfectly dithered frames, the performance of drizzle is spectacular in this test, being able to recover the data destroyed by sparse sampling at 2:1 scale. Don't let this result convince you that Bayer drizzle x2 is so easy in practice, however. This is a synthetic test, not real data, dithering is ideal, and we are speaking of 100 frames.

I have made many tests like this one with sets of 100 randomly dithered frames, with nearly identical results, as expected. For some extremely weird reason that I can't understand, the result after fixing the bug is worse with the test set that I generated for my initial post:

However, it is fantastic with drizzle x2, as you can see above. The only thing that changes among different data sets is the fractional dithering displacements. A set of 100 frames is not a large sample in statistical terms, and chances are that a particular combination of displacements are suboptimal for drizzle. It is not that the drizzle result is now a bad one---it is quite good, actually, with minimal color artifacts and all pattern structures reconstructed quite well---, but the result without fixing the bug was so good...

By the way, I have now a new version of the synthetic data generation script, which I am using to perform intensive testing of different tools and new pixel rejection algorithms:

Code: [Select]

/*
 * A synthetic data generation script. Useful for testing image integration
 * and rejection algorithms.
 *
 * Written by Juan Conejero, PTeam.
 * October 2016
 */

// View identifier of the original image.
#define ORIGINAL_IMAGE_ID     "original"

// Number of synthetic images to generate.
#define NUMBER_OF_IMAGES      100

// Output base directory.
#define OUTPUT_BASE_DIR       "/home/juan/tmp/bayer-drizzle-test/3"

// Generate synthetic images dithered by integer pixel amounts?
#define INTEGER_DITHERING     false

// Generate synthetic images dithered by fractional pixel amounts?
#define FRACTIONAL_DITHERING  true

// Maximum image displacement in pixels (absolute value).
#define MAX_DELTA_PX          25

// Probability of occurrence of cosmic rays in an image.
#define PROB_COSMICS          0.25

// Maximum number of artificial cosmic rays in an image.
#define MAX_COSMICS           4

// Minimum fraction of artificial outlier pixels.
#define MIN_OUTLIERS          0.0001

// Maximum fraction of artificial outlier pixels.
#define MAX_OUTLIERS          0.001

// Generate Bayer mosaiced data?
#define MOSAICED_DATA         true

// Generate Bayer RGB or monochrome CFA images?
#define BAYER_RGB             false

// Bayer pattern for sparse sampling of synthetic data, when MOSAICED_DATA is
// true. Can be one of "RGGB", "BGGR", "GRBG", and "GBRG".
#define BAYER_PATTERN         "RGGB"

#include <pjsr/Interpolation.jsh>
#include <pjsr/UndoFlag.jsh>

function translateImageByIntegerPixels( image, dx, dy )
{
   image.shiftBy( Math.round( dx ), Math.round( dy ) );
}

function translateImageByFractionalPixels( image, dx, dy )
{
   image.translate( -dx, -dy );
}

function generateOutlierPixels( image, fraction )
{
   let count = Math.round( fraction * image.numberOfPixels );
   for ( let z = 0; z < image.numberOfChannels; ++z )
      for ( let i = 0; i < count; ++i )
      {
         let x = Math.trunc( image.width*Math.random() );
         let y = Math.trunc( image.height*Math.random() );
         image.setSample( Math.random(), x, y, z );
      }
}

function generateCosmics( image, count )
{
   if ( count > 0 )
   {
      let b = new Bitmap( image.width, image.height );
      b.fill( 0x00000000 );
      let g = new VectorGraphics( b );
      g.antialiasing = true;
      for ( let i = 0; i < count; ++i )
      {
         let x1 = Math.trunc( image.width*Math.random() );
         let y1 = Math.trunc( image.height*Math.random() );
         let x2 = Math.trunc( image.width*Math.random() );
         let y2 = Math.trunc( image.height*Math.random() );
         g.pen = new Pen( 0xffffffff, 1 + 2*Math.random() );
         g.drawLine( x1, y1, x2, y2 );
      }
      g.end();

      image.blend( b );
   }
}

function bayerChannelCoordinates( bayerPattern )
{
   let Rx, Ry, Gx0, Gx1, Bx, By;
   switch ( bayerPattern )
   {
   // R G
   // G B
   default:
   case "RGGB":
      Rx  = 0; Ry  = 0;
      Gx0 = 1; Gx1 = 0;
      Bx  = 1; By  = 1;
      break;
   // B G
   // G R
   case "BGGR":
      Rx  = 1; Ry  = 1;
      Gx0 = 1; Gx1 = 0;
      Bx  = 0; By  = 0;
      break;
   // G R
   // B G
   case "GRBG":
      Rx  = 1; Ry  = 0;
      Gx0 = 0; Gx1 = 1;
      Bx  = 0; By  = 1;
      break;
   // G B
   // R G
   case "GBRG":
      Rx  = 0; Ry  = 1;
      Gx0 = 0; Gx1 = 1;
      Bx  = 1; By  = 0;
      break;
   }
   return { Rx:Rx, Ry:Ry, Gx0:Gx0, Gx1:Gx1, Bx:Bx, By:By };
}

function convertRGBImageToBayerRGB( view, bayerPattern )
{
   let c = bayerChannelCoordinates( bayerPattern );
   let P = new PixelMath;
   P.expression0 = "iif( x()%2 == " + c.Rx.toString() + " && y()%2 == " + c.Ry.toString() + ", $T, 0 )";
   P.expression1 = "iif( iif( y()%2 == 0, x()%2 == " + c.Gx0.toString() + ", x()%2 == " + c.Gx1.toString() + " ), $T, 0 )";
   P.expression2 = "iif( x()%2 == " + c.Bx.toString() + " && y()%2 == " + c.By.toString() + ", $T, 0 )";
   P.useSingleExpression = false;
   P.createNewImage = false;
   P.executeOn( view, false/*swapFile*/ );
}

function convertRGBImageToBayerCFA( view, bayerPattern )
{
   let c = bayerChannelCoordinates( bayerPattern );
   let P = new PixelMath;
   P.expression0 = "iif( y()%2 == " + c.Ry.toString() + ", iif( x()%2 == " + c.Rx.toString() + ", $T[0], $T[1] ), iif( x()%2 == " + c.Bx.toString() + ", $T[2], $T[1] ) )";
   P.expression1 = "0";
   P.expression2 = "0";
   P.useSingleExpression = false;
   P.createNewImage = false;
   P.executeOn( view, false/*swapFile*/ );
   view.beginProcess( UndoFlag_NoSwapFile );
   let image = Image.newUIntImage( 16 );
   view.image.selectedChannel = 0;
   image.assign( view.image );
   view.image.resetChannelSelection();
   view.image.assign( image );
   view.endProcess();
}

function main()
{
   if ( !File.directoryExists( OUTPUT_BASE_DIR ) )
      File.createDirectory( OUTPUT_BASE_DIR );

   let outputIntegerDir = OUTPUT_BASE_DIR + "/integer";
   let outputIntegerBayerDir = OUTPUT_BASE_DIR + "/integer/bayer";
   if ( INTEGER_DITHERING )
   {
      if ( !File.directoryExists( outputIntegerDir ) )
         File.createDirectory( outputIntegerDir );
      if ( MOSAICED_DATA )
         if ( !File.directoryExists( outputIntegerBayerDir ) )
            File.createDirectory( outputIntegerBayerDir );
   }

   let outputFractionalDir = OUTPUT_BASE_DIR + "/fractional";
   let outputFractionalBayerDir = OUTPUT_BASE_DIR + "/fractional/bayer";
   if ( FRACTIONAL_DITHERING )
   {
      if ( !File.directoryExists( outputFractionalDir ) )
         File.createDirectory( outputFractionalDir );
      if ( MOSAICED_DATA )
         if ( !File.directoryExists( outputFractionalBayerDir ) )
            File.createDirectory( outputFractionalBayerDir );
   }

   let originalWindow = ImageWindow.windowById( ORIGINAL_IMAGE_ID );

   for ( let i = 0; i < NUMBER_OF_IMAGES; ++i )
   {
      let window = new ImageWindow( 1, 1, 1, 16/*bitsPerSample*/, false/*floatSample*/ );
      let view = window.mainView;

      let dx = Math.trunc( MAX_DELTA_PX*Math.random() ) + Math.random();
      if ( Math.random() < 0.5 )
         dx = -dx;

      let dy = Math.trunc( MAX_DELTA_PX*Math.random() ) + Math.random();
      if ( Math.random() < 0.5 )
         dy = -dy;

      let fractionOfOutliers = MIN_OUTLIERS + (MAX_OUTLIERS - MIN_OUTLIERS)*Math.random();

      let numberOfCosmics = (Math.random() < PROB_COSMICS) ? Math.round( MAX_COSMICS*Math.random() ) : 0;

      if ( INTEGER_DITHERING )
      {
         view.beginProcess( UndoFlag_NoSwapFile );
         view.image.assign( originalWindow.mainView.image );
         translateImageByIntegerPixels( view.image, dx, dy );
         generateOutlierPixels( view.image, fractionOfOutliers );
         generateCosmics( view.image, numberOfCosmics );
         view.endProcess();

         if ( MOSAICED_DATA )
            if ( BAYER_RGB )
               convertRGBImageToBayerRGB( view, BAYER_PATTERN );
            else
               convertRGBImageToBayerCFA( view, BAYER_PATTERN );

         window.saveAs( (MOSAICED_DATA ? outputIntegerBayerDir : outputIntegerDir) + format( "/test%04d.xisf", i+1 ),
                        false/*queryOptions*/, false/*allowMessages*/, true/*strict*/, false/*verifyOverwrite*/ );
      }

      if ( FRACTIONAL_DITHERING )
      {
         view.beginProcess( UndoFlag_NoSwapFile );
         view.image.assign( originalWindow.mainView.image );
         view.image.interpolation = Interpolation_Bilinear;
         translateImageByFractionalPixels( view.image, dx, dy );
         generateOutlierPixels( view.image, fractionOfOutliers );
         generateCosmics( view.image, numberOfCosmics );
         view.endProcess();

         if ( MOSAICED_DATA )
            if ( BAYER_RGB )
               convertRGBImageToBayerRGB( view, BAYER_PATTERN );
            else
               convertRGBImageToBayerCFA( view, BAYER_PATTERN );

         window.saveAs( (MOSAICED_DATA ? outputFractionalBayerDir : outputFractionalDir) + format( "/test%04d.xisf", i+1 ),
                        false/*queryOptions*/, false/*allowMessages*/, true/*strict*/, false/*verifyOverwrite*/ );
      }

      window.forceClose();

      console.noteln( format( "* %4d : dx=%+.2lf dy=%+.2lf o=%.4lf c=%d", i+1, dx, dy, fractionOfOutliers, numberOfCosmics ) );
   }
}

main();

The script is now much more flexible. It generates monochrome Bayer CFA frames by default, and can generate artificial outliers and simulated cosmic rays. BTW, this script has allowed me to discover very small pixel rejection errors in the DrizzleIntegration tool. This is a different (very minor) bug that I'll comment on in a different post.

However I completely understand if you have already spent enough time experimenting with Bayer Drizzle

I cannot get tired of experimenting and finding and fixing bugs and issues; I simply cannot afford that

Do you have a GitHub account, or can create one (it's free)? I'd be glad to have you in our GitHub team. Just let me know your GitHub user name and I'll invite you, if you wish. I make this extensive to anyone interested in PixInsight development.

[sharkmelley]

Hi Juan,

Thanks for the comprehensive reply and for your effort on this. I'm certainly glad there will be a bug fix on the way.

Yes, I'm definitely interested in PixInsight development and I already have a few ideas in prototype form as PI scripts.  I'll send you a GitHub user name.  I did download the zip last night and built it (MS Visual Studio Express 2013 on W7) but was unable to install the resulting ImageIntegration-pxm.dll into PixInsight.  I wonder if it is because the post-build "signtool" step failed?  I'll try again and start a separate thread on PCL and PJSR Development if I still have problems.

[Later comment]  Above problem now solved - I mistakenly compiled as Win32 instead of x64.  With x64 I can now install it in PI and debug it.

Mark

[sharkmelley]

Hi Juan,

I have built a version of ImageIntegration with your fixes and it is working really well.  However, I have a couple of observations for you. 

Firstly, I think I've found the reason for the earlier puzzling runs, where some results were better than others.  It all depends on which registration reference image is selected.  If the reference image is approximately an integer offset from the original image then the (bug fixed) drizzle does an excellent job of reproducing all the features in the original test image.  However, if the reference image is offset by approximately 0.5 pixels then drizzle is registering to something already interpolated and this will smear the features of the stacked image when compared to the original test image.  These are the two extremes - in general the reference image will be somewhere between these extremes.  I think your data generation script needs a option to generate "image1" with zero offset.  This can then be used as an optimal registration reference image for testing purposes.

The second observation concerns the drizzle files.  When the BatchPreprocessing script is used with "Generate drizzle data" switched on, it produces a set of registered images and drizzle files in the "registered" folder and a set of bayer drizzle files in the "registered/bayer" folder.  Now if I subsequently run ImageIntegration, I "Add Files" from the "registered" folder but then "Add Drizzle Files" automatically displays the drizzle files from the same "registered" folder. It is very easy to mistakenly select these but they are the wrong ones to use for bayer drizzle.  For bayer drizzle I need to pick up the drizzle files from the "registered/bayer" folder instead.  This has caught me out on numerous occasions and it gives an inferior result i.e. not bayer drizzled.  I wonder how users can be helped to avoid this trap?

Mark