Some one asked me how I choose the matrix coef to convert a five multispectral image to a three (RGB) multispectral image. I feel that there is a best way to do that by maximising some distance in a color space. How would you do that? Examples of multispectral images are available here
More details on the input images would be helpful: are they the result of performing spectral unmixing, so each image represents the abundance of a particular fluorophore? Or are the images the raw results? The best approach would probably differ between the two options. If the images are unmixed, then you can calculate the "true" spectrum at each pixel using the images and the emission spectra from the fluorophores. Then using these spectra, you could calculate the effective color in a number of ways, probably the simplest being a weighted average of the "red", "green", and "blue" regions of the spectra. You'd want to look at the psychophysics literature to figure out how best to do this, but I'm sure there are a ton of references online for converting a spectrum into an RGB color. If the images are mixed (that is, raw), it's even easier: figure out the "apparent color" of the emission filter for each image (see below for sample code), and then colorize each image by the apparent color, and screen the images together. (Code for screening a set of images together is below, too -- 'multi_screen()') Or, returning to the unmixed case, you could forget using the "natural" colors, and instead choose a set of colors maximally dissimilar from one another in some color space as the color for each fluorophore. This might make visualization simpler. You'd choose the colors (not sure how best), and then colorize the images and screen them together. Zach def screen(src, dst, max=255): return max - (((max - src)*(max - dst))/max) def multi_screen(images, max=255): dst = images[0] for image in images[1:]: dst = screen(image, dst, max) return dst def wavelength_to_rgb(l): # http://www.physics.sfasu.edu/astro/color/spectra.html # So-called "Bruton's Algorithm" assert (350 <= l <= 780) if l < 440: R = (440-l)/(440-350.) G = 0 B = 1 elif l < 490: R = 0 G = (l-440)/(490-440.) B = 1 elif l < 510: R = 0 G = 1 B = (510-l)/(510-490.) elif l < 580: R = (l-510)/(580-510.) G = 1 B = 0 elif l < 645: R = 1 G = (645-l)/(645-580.) B = 0 else: R = 1 G = 0 B = 0 if l > 700: intensity = 0.3 + 0.7 * (780-l)/(780-700.) elif l < 420: intensity = 0.3 + 0.7 * (l-350)/(420-350.) else: intensity = 1 return intensity*numpy.array([R,G,B], dtype=float)