Alternative Debayer algorithm

fredvanner

Well-known member
The supported algorithms all have issues for astronomical use:
bilinear: interpolation smooths raw data, suppressing local detail.
VNG: fine for photography because preserves discontinuous boundaries, but not good for astronomy - stars are surrounded by a discontiuity, so VNG smooths across the interior of stars, flattening their profile (and thus modifying their PSF - e.g. for convolution).
superpixel: discards the higher resolution available in the green channel data, and , worse, averages the green channel, losing the local luminosity maxima that contribute to detail (I think max(G1, G2) would be better for astronomy).
All three methods end up smoothing the valuable, higher resolution green channel luminosity data.
I have tried an alternative artefact-free method. It is basically a sliding superpixel window, but selecting only one of the two possible green values. Thus each pixel is constructed from exactly one R, one G and one B from the original RGGB image. Each R and B is used in four separate RGB pixels; each G is used in two RGB pixels. This process has the advantage that every RGB value corresponds to an original image value (no averaging / smoothing / interpolation, and all RGB values, including maxima, are preserved), and that that the higher resolution of the green channel is not lost in the RGB image.

I would attach an example, but the files are too large (and compressed versions, such as jpg, would defeat the object); however I attach a jpg of a sample debayered image (M51 240s ASI183MC Pro, dark subtracted).

Curently, I have implemented this with PixelMath, but as a complete PixInsight beginner, I can't work out how to put it together as a process that I can apply to lists of images (basic ProcessContainer approaches stall because I can't see how to create a new blank RGB view for the RGB merge; more flexible scripting approaches stall because I can find no documentation for the js API; any help / pointers to documentation would be much appreciated.)
 

Attachments

  • Test_newRGB.jpg
    497.8 KB · Views: 119
I believe this algorithm is better for use by astrophotographers (while entirely inappropriate for terrestrial photography). As an artefact-free, loss-free (no original data values are modified) algorithm, it leaves the decision of where in the image processing sequence to perform any smoothing. In contrast, biliner and VNG processing smooths the raw data, removing local maxima and minima. Not surprisingly, the differences are small and subtle, but I believe they are in the spirit of most astrophotographers: "I want to decide what happens to my raw data!".
I would really appreciate a second opinion, particularly from a more experienced astrophotographer. If anyone would like to point me at a raw RGGB file, I would be happy to publish my debayered RGB version.
(BTW, I'm not a "lunatic fringe" guy. I'm a professional engineer with a maths degree, so I do understand what I'm doing ... I think!)
 
If you want to use a JS script then you have to look at existing scripts and work it out, sorry but this is the only way at the moment, it is the way we have all done it.....plus the odd prompts from Juan when we get stuck.
 
Fred, please post your pixelmath if possible. I might be able to write a very basic starting script that you would be able to work with.
 
Fred, here is an attempt. Of course you can use PI's script editor to create, edit, and execute this .js file. You can run it either on the active window or a set of files (see end of code). The arrays nextPattern and delta are central, to phase across the raw and to choose components properly. My choices here are likely not what you desire. But maybe you can use this as a starting point and make it work.

JavaScript:
// AlternativeDebayer.js

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

// Gives a new image window.
function newImageWindow(image, fullId) {
   var imageWindow = new ImageWindow(
      image.width,
      image.height,
      image.numberOfChannels,
      image.bitsPerSample,
      image.isReal,
      image.numberOfChannels > 1,
      fullId
   );

   imageWindow.mainView.beginProcess(UndoFlag_NoSwapFile);
   imageWindow.mainView.image.assign(image);
   imageWindow.mainView.endProcess();

   return imageWindow;
};

// Imports an image from filePath.
function importImage(filePath) {
   var imageWindows = ImageWindow.open(filePath);
   if (imageWindows.length == 0) {
      throw "missing file";
   }

   var image = new Image(imageWindows[0].mainView.image);

   for (var i = 0; i != imageWindows.length; ++i) {
      imageWindows[i].forceClose();
   }

   return image;
}

// Exports an image to filePath.
function exportImage(filePath, image, overwrite) {
   var imageWindow = newImageWindow(image, "export");
   imageWindow.saveAs(filePath, false, true, true, !overwrite);
   imageWindow.forceClose();

   return image;
}

// Bayer patterns
function Bayer() {
}
Bayer.RGGB = 0;
Bayer.BGGR = 1;
Bayer.GBRG = 2;
Bayer.GRBG = 3;

function AlternativeDebayer() {
   // Gives next pattern for each phase.
   // invariant nextPattern[p][0][0] == p for all pattern p.
   this.nextPattern = [
      [
         [Bayer.RGGB, Bayer.GRBG],
         [Bayer.GBRG, Bayer.BGGR]
      ],
      [
         [Bayer.BGGR, Bayer.GBRG],
         [Bayer.GRBG, Bayer.RGGB]
      ],
      [
         [Bayer.GBRG, Bayer.BGGR],
         [Bayer.RGGB, Bayer.GRBG]
      ],
      [
         [Bayer.GRBG, Bayer.RGGB],
         [Bayer.BGGR, Bayer.GBRG]
      ]
   ];

   // Gives delta [dx, dy] for each channel. May be -1, 0, or 1.
   this.delta = [
      [
         [0, 0], [0, 1], [1, 1]
      ],
      [
         [1, 1], [0, 1], [0, 0]
      ],
      [
         [0, 1], [0, 0], [1, 0]
      ],
      [
         [1, 0], [0, 0], [0, 1]
      ]
   ];

   // Gives reflection of delta d along edge of image.
   this.reflectDelta = function(d, z, size) {
      return z == 0 ? Math.abs(d) : z == size - 1 ? -Math.abs(d) : d;
   };

   // Gives debayer of raw wrt pattern.
   this.debayerRaw = function(raw, pattern) {
      var rgb = new Image(
         raw.width,
         raw.height,
         3,
         ColorSpace_RGB,
         raw.bitsPerSample,
         raw.sampleType
      );

      var dx;
      var dy;
      var v = 0;
      var p = pattern;
      for (var y = 0; y != raw.height; ++y) {
         var u = 0;
         p = this.nextPattern[pattern][v][u];
         for (var x = 0; x != raw.width; ++x) {
            dx = this.reflectDelta(this.delta[p][0][0], x, raw.width);
            dy = this.reflectDelta(this.delta[p][0][1], y, raw.height);
            rgb.setSample(raw.sample(x + dx, y + dy, 0), x, y, 0);
            dx = this.reflectDelta(this.delta[p][1][0], x, raw.width);
            dy = this.reflectDelta(this.delta[p][1][1], y, raw.height);
            rgb.setSample(raw.sample(x + dx, y + dy, 0), x, y, 1);
            dx = this.reflectDelta(this.delta[p][2][0], x, raw.width);
            dy = this.reflectDelta(this.delta[p][2][1], y, raw.height);
            rgb.setSample(raw.sample(x + dx, y + dy, 0), x, y, 2);

            u = 1 - u;
            p = this.nextPattern[pattern][v][u];
         }

         v = 1 - v;
         p = this.nextPattern[pattern][v][u];
      }

      return rgb;
   };
}

// Debayer the active window wrt pattern.
function debayerActiveWindow(pattern, suffix) {
   var mainView = ImageWindow.activeWindow.mainView;
   if (!mainView.isView) {
      throw "no active window";
   }
   var fullId = mainView.fullId;

   var raw = new Image(mainView.image);
   if (raw.numberOfChannels != 1) {
      throw "not monochannel raw";
   }
   if (raw.width < 2 || raw.height < 2) {
      throw "raw too small";
   }

   var alternativeDebayer = new AlternativeDebayer();
   var rgb = alternativeDebayer.debayerRaw(raw, pattern);

   var imageWindow = newImageWindow(rgb, fullId + suffix);
   imageWindow.show();

   raw.free();
   rgb.free();
}

// Debayer an image from filePath wrt pattern.
function debayerFilePath(filePath, pattern, suffix, overwrite) {
   var raw = importImage(filePath);
   if (raw.numberOfChannels != 1) {
      throw "not monochannel raw";
   }
   if (raw.width < 2 || raw.height < 2) {
      throw "raw too small";
   }

   var alternativeDebayer = new AlternativeDebayer();
   var rgb = alternativeDebayer.debayerRaw(raw, pattern);

   exportImage(File.appendToName(filePath, suffix), rgb, overwrite);

   raw.free();
   rgb.free();
}

// Debayer the active window wrt pattern.
if (true) {
   var pattern = Bayer.RGGB;
   var suffix = "_d";
   debayerActiveWindow(pattern, suffix);
}

// Debayer a set of images from filePaths wrt pattern.
if (false) {
   var pattern = Bayer.RGGB;
   var suffix = "_d";
   var overwrite = true;
   var directoryPath = "/Users/Home/Documents/Astronomy/PixInsight/AlternativeDebayer/";
   var filePaths = [
      "light_01.fit",
      "light_02.fit",
      "light_03.fit"
   ];
   for (var i = 0; i != filePaths.length; ++i) {
      debayerFilePath(directoryPath + filePaths[i], pattern, suffix, overwrite);
   }
}
 
Thanks a lot. I'm a beginner withPixinsight and js, but I've spent 50 years programming in everything from FORTRAN to C++, so I'll work it out given time! As requested, I attach the basic PixelMath code for the algorithm (I haven't yet bothered to tidy up the last row and column, where it falls off the edge - but this is no big issue). Oops - it won't let me attach a .js file; how do I do that?
 
OK, found "insert code" under the "..." dropdown!
I will try and sort it out myself, but for the record the basic PixelMath (probably easily improved) is:
Code:
var P = new PixelMath;
P.expression = "P00=pixel($T,x(), y()); P01=pixel($T,x(),y()+1); P10=pixel($T, x()+1, y()); P11=pixel($T, x()+1, y()+1); x0=frac(x()/2)==0; y0=frac(y()/2)==0; iif(x0,iif(y0, P00, P01),iif(y0, P10, P11))";
P.expression1 = "P00=pixel($T,x(), y()); P01=pixel($T,x(),y()+1); P10=pixel($T, x()+1, y()); P11=pixel($T, x()+1, y()+1); x0=frac(x()/2)==0; y0=frac(y()/2)==0; iif(x0,iif(y0, P10, P00),iif(y0, P00, P10))";
P.expression2 = "P00=pixel($T,x(), y()); P01=pixel($T,x(),y()+1); P10=pixel($T, x()+1, y()); P11=pixel($T, x()+1, y()+1); x0=frac(x()/2)==0; y0=frac(y()/2)==0; iif(x0,iif(y0, P11, P10),iif(y0, P01, P00))";
P.expression3 = "";
P.useSingleExpression = false;
P.symbols = "P00, P01, P10, P11, x0, y0";
P.generateOutput = true;
P.singleThreaded = false;
P.use64BitWorkingImage = false;
P.rescale = false;
P.rescaleLower = 0;
P.rescaleUpper = 1;
P.truncate = true;
P.truncateLower = 0;
P.truncateUpper = 1;
P.createNewImage = true;
P.showNewImage = true;
P.newImageId = "";
P.newImageWidth = 0;
P.newImageHeight = 0;
P.newImageAlpha = false;
P.newImageColorSpace = PixelMath.prototype.RGB;
P.newImageSampleFormat = PixelMath.prototype.f32;
 
Ok, on comparison, results seem same on a test natural image, except for anisotropy axis. Your's has more zippering on one axis, mine on the other axis, i.e., the same up to a symmetry. This is due to differing pixel choice coded in the delta table, I expected differences like this to be likely.

Clearly, you could code in your expressions and remove the table lookup, which would likely make things easier for you, but you would of course loose the automatic edge handling and support for the other bayer patterns, which you might want to add, probably easily.
 
Last edited:
Yup. The 1x2 G array can be vertical or horzontal with no preference (you have to break the symmetry). I think I follow your code (including misplaced closing brace in the Bayer() function), but theres lots of new Pixinsight stuff (understanding the available objects, their properties and methods) that is difficult to extract for scripts; some basic documentation would really help - like what global objects are available to a script without declaration.
 
Re: symmetry, yes, with the delta table below (flipped choice) the results are identical.

JavaScript:
this.delta = [
      [
         [0, 0], [1, 0], [1, 1]
      ],
      [
         [1, 1], [1, 0], [0, 0]
      ],
      [
         [0, 1], [0, 0], [1, 0]
      ],
      [
         [1, 0], [0, 0], [0, 1]
      ]
   ];

Regarding documentation, have you looked at PI's Object Explorer? Also, .jsh headers in your PI installation folder include/pjsr.
 
Last edited:
Thanks for your time and effort.
Yes, I had looked at the object browser, but the list of properties and methods with no further information (what do the methods do? how do you call them? what are the property types? what do the properties contain?) is simply frustrating. I had looked at some of the headers in include/pjsr, but they had mostly contained definitions for the names of integer-coded arguments - which aren't much use if you don't know the names in the first place! However, on further inspection there are some .jsh files which do contain the calling interface of some object methods. It looks as though the only effective way ahead, as you first suggested, is to "reverse engineer" enough of the supplied scripts. That sounds like a project for several months, and I'm not sure I have the stamina! ... but maybe during COVID lockdown I will find the time!
I've tried your script and it works fine - and resolves (by reflection) the boundary problems that I hadn't bothered to fix. I think that will be my preferred DB method in future (that way I can do the smoothing when and how I want).
 
Glad it works. Sorry for the lack of documentation. This is a problem.

Rather than editing file paths as strings in the code, the code below makes the script easier to use with an OpenFileDialog to let you select files for debayer.

Of course, you can use this entire code as an example template for other processing operations.

JavaScript:
// Debayer a set of images from filePaths wrt pattern.
if (true) {
   var pattern = Bayer.RGGB;
   var suffix = "_d";
   var overwrite = true;

   var openFileDialog = new OpenFileDialog;
   openFileDialog.multipleSelections = true;
   openFileDialog.caption = "Select frames for debayer";
   openFileDialog.filters = [
      ["All supported formats", ".fit", ".fits", ".fts", ".xisf"],
      ["FITS Files", ".fit", ".fits", ".fts"],
      ["XISF Files", ".xisf"]
   ];
   if (openFileDialog.execute()) {
      for (var i = 0; i != openFileDialog.fileNames.length; ++i) {
         debayerFilePath(openFileDialog.fileNames[i], pattern, suffix, overwrite);
      }
   }
}
 
Thanks! I had just spent a day trying to do this file based option (by trying to disentangle the file interface bits from the BCE script). It doesn't help that I'm also a beginner with js (I'm still trying to get my head round the syntax of:
JavaScript:
// Bayer patterns
function Bayer() {
}
Bayer.RGGB = 0;
Bayer.BGGR = 1;
Bayer.GBRG = 2;
Bayer.GRBG = 3;
in your script - it seems there would be lots of more straightforward ways to achieve the same result, but I'm probably missing some hard-learned subtlety).
Thanks again for your time.
 
Sorry for the confusion.

I envision static properties so that their value is available without instancing, in a sort of base class with nothing else defined yet.
 
Hi Fred, Mike,

Is an implementation of this algorithm available through the PixInsight update mechanism?

Thanks,
John
 
I've learned a bit more about astro processing (and PI) since I originally proposed this. I have slightly improved (I believe) the algorithm, getting rid of the highly visible asymmetry (effectively by alternating the "horizontal" and "vertical" G channel allocation between 2x2 Bayer cells). I attach a PixelMath implementation. I have not pushed ahead with a script version since I now understand that a proper "Debayer" implementation needs to update the CFA path property (so subsequent StarAlignment can build the drizzle file if requested). I haven't yet researched this enough to implement it with confidence.
I still believe that this method is much better than VNG for astrophotography. In most cases it will be very similar to "bilinear", but definitely superior to "superpixel" (since it retains the higher luminosity resolution, particularly in the G channel).
 

Attachments

  • nnDebayer.zip
    999 bytes · Views: 80
... further comment: It was only after further reading that I realised that my prosed method is effectively the same as the "nearest neighbour" debayer algorithm - the grandfather of all debayer algorithms. For terrestrial photography a number of more "advanced" algorithms were developed, mainly to protect the sharp, locally linear, gradient discontinuities that occur in many "artificial" environments (e.g. rectilinear architecture). For astrophotography, where the main steep local gradients are associated with stars, and are radial rather than linear in structure, these algorithms have a strong tendency to produce uninteded (if small and very localised) artefacts. This is best avoided by unbiased (i.e. bilinear) interpolation, or no interpolation at all ("nearest neighbour").
 
There is one subtle defect (feature) of this method that I have discovered. Because the "B" cell of the bayer matrix is used as the "B" value for the whole 2x2 cell, it's centre is effectively shifted up and left by half a pixel (assuming RGGB); similarly, the "R" cell is effectively shifted down and right by half a pixel. This can be compensated by using ChannelMatch to shift R by (-0.7,-0.7) and B by (0.7,0.7) (0.7 should actually be 0.7071... or 1/root(2)) (-0.5,-0.5) and B by (0.5,0.5) - the distance shifted is 0.7071, but the (x,y) shift is +/- (0.5,0.5) - my mistake. It is very hard to see the difference!
Note that this is not interpolation between separate "B" or "R" values, it is a bulk shift of the whole "B" and "R" layer.
The following shows the difference:
Left: bilinear; Middle: nearest neighbour - no shift; Right: nearest neighbour with R & B compensating shift.
1616263357505.png

Note the lack of interpolation (smoothing) with nn. Why not let your integration smooth the more precise nn data.
Any implementation of this algorithm needs to include this compensation (as well as handling all the other bayer matrix combinations).
 
Last edited:
Back
Top