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 '14 at 9:33
Why not ColorCombine[Image/@({m13R, m13G, m13B}-dark-bias), "RGB"]? –  Kuba May 1 '14 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. –  nbogs May 1 '14 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. –  Simon Woods May 1 '14 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 –  nbogs May 1 '14 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"][[1]] & /@ Rest@FileNames["*CCD*.fit"]; (* 3 pics *)

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

dark = Import[FileNames["*dark*.fit"][[1]], "RawData"][[1]]; (* 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. –  Oleksandr R. May 2 '14 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 '14 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. –  Wizard May 1 '14 at 11:37