# Combining three R, G, B matrices into a single matrix of 3-tuples?

Working with Mathematica 8 on OSX here.

I have three $1024 \times 1024$ matrices representing the counts from consecutive exposures onto a CCD with $R$, $G$, and $B$ filters on it, and am looking to combine the three matrices to get a nice color picture.

My first shot was something of the following sort:

ColorCombine[{
MatrixPlot[m13R - dark - bias, ColorFunction -> (RGBColor[#, 0, 0] &)],
MatrixPlot[m13G - dark - bias, ColorFunction -> (RGBColor[0, #, 0] &)],
MatrixPlot[m13B - dark - bias, ColorFunction -> (RGBColor[0, 0, #] &)]
}]


but it looks like ColorCombine messes with the resolution, and so what was originally $1024\times1024$ comes out looking rather pixelated.

What I'm thinking now is to do the equivalent of Transpose but am unsure as to how to combine my three matrices into a $1024\times1024$ matrix with a $3$-tuple at each entry. From there I think Image[combined, ColorSpace -> "RGB"] will work to plot the now-colored image.

For those interested, you can find the relevant data here.

Any help appreciated!

• You can use MapThread[List, {m13R, m13G, m13B}, 2] (assuming these are n x m matrices) and then ArrayPlot on the output.
– gpap
May 1, 2014 at 9:33
• Why not ColorCombine[Image/@({m13R, m13G, m13B}-dark-bias), "RGB"]?
– Kuba
May 1, 2014 at 9:43
• @gpap, using MapThread[List, {m13R - k, m13G - k, m13B - k}, 2], where $k=dark + bias$, in conjunction with Image[-2724/RGB, ColorSpace -> "RGB"] provides a somewhat reasonable picture--at least one that can be adjusted to look good. ($-2724$ is a normalizing factor to scale matrix values into the range $0-1$.) @kuba, @SimonWoods: Unfortunately both Image[{m13R/65492, m13G/65530, m13B/65518}, Interleaving -> False] and ColorCombine[Image/@({m13R, m13G, m13B}-dark-bias), "RGB"] produce black either black frames or a list of images, respectively. May 1, 2014 at 18:30
• If you are getting a black frame it suggests the arrays are scaled incorrectly. Try using Rescale[{r, g, b}] or using ImageAdjust on the final image. May 1, 2014 at 18:58
• @Kuba I just checked the filter used on "Copy of CCD Image 67.fit" and for me it shows up red, as expected. The filter comes up upon doing Import[path, "Metadata"] // TableForm May 1, 2014 at 21:30

SetDirectory@NotebookDirectory[];
FileNames["*.fit"]

{"Copy of bias6.fit", "Copy of bias7.fit", "Copy of bias8.fit", "Copy of bias9.fit",
"Copy of CCD Image 64.fit", "Copy of CCD Image 65.fit", "Copy of CCD Image 66.fit",
"Copy of CCD Image 67.fit", "Copy of dark12.fit"}


pics = Import[#, "RawData"][] & /@ Rest@FileNames["*CCD*.fit"]; (* 3 pics *)

bias = Import[#, "RawData"][] & /@ FileNames["*bias*.fit"]; (*4 pics *)

dark = Import[FileNames["*dark*.fit"][], "RawData"][]; (* 1 pic *)


### Data info:

{#, "EXPTIME", "FILTER"} /. Import[#, "Metadata"] & /@ FileNames["*.fit"] // Grid ### Processing

Exposure time is important, we have to rescale dark to the exposure time of observations. Here it means we have to divide by 60/20 == 3.

Notice that dark has the bias so we have to subtract it before rescaling to not rescale bias. We can subtract bias from observations separately later.

reducedDark = (dark - Mean[bias]) /3. ;

reducedObs = (# - Mean@bias - reducedDark) & /@ pics;
{mean, std} = {Mean@#, StandardDeviation@#} &@
DeleteCases[ Clip[Flatten@reducedObs, {280, ∞}, {0, ∞}], 0]

{331.636, 402.461}


And using, Simon's suggestion:

ImageAdjust[Image[Reverse@reducedObs, Interleaving -> False],
{0, 0, 1}, mean + {-.3 std, 1.5 std}, {0, 1}] Of course it's not even close to the observations reduction in astronomy. What we are missing here are e.g.:

etc. etc. I'm not pretending I know what they do :) I've just faced an idea of this on lab couple years ago. :)

• Note that to get the colors correct also requires consideration of the relative sensitivity in each band. Incidentally, RawData can often be integer ADC counts, in which case arithmetic manipulation of the images might be painfully slow. This didn't occur in version 8 only because it had a bug whereby the data were returned as floating-point values even if the file only contained integers. May 2, 2014 at 20:11
• @OleksandrR. You mean CCD sensitivity right? Also filter profiles are not just top hats what readers should keep in mind.
– Kuba
May 3, 2014 at 10:19

You can use Transpose with its second argument:

(* example data *)
{r, g, b} = DiskMatrix[#, 100] & /@ {30, 20, 10};

Image[Transpose[{r, g, b}, {3, 1, 2}]] But it is simpler to use the Interleaving option of Image:

Image[{r, g, b}, Interleaving -> False] • Didn't know about the Interleaving Option yet. Good to know. May 1, 2014 at 11:37