Thursday, February 25, 2010

Transforming RGB data to wavelength information

Here, I'm presenting a reverse transformation to wavelength to RGB. Right now I have no means of testing its accuracy on some real data. As a first check, I've run it on the mercury lamp image shown in my older post. There are small deviations caused probably by low intensity of light in some places. If used on some real stuff, I would definitely attempt to first fit a line (or some other meaningful function) through the estimated wavelengths. Then of course plot the spectra against fitted values. The domain of this function is limited to 380 - 645 nm. So, here's the algorithm in Python:

def RGB2lambda(R, G, B):
"""Returns 0 if indeciferable"""
# selects range by maximum component
# if max is blue - range is 380 - 489
# if max is green - range is 490 - 579
# if max is red - range is 580 - 645

# which colour has highest intensity?
high = float(R)
highind = 1
if G > high:
high = float(G)
highind = 2
if B > high:
high = float(B)
highind = 3

# normalize highest to 1.0
RGBnorm = [R / high, G / high, B / high]

# start decoding
RGBlambda = 0
if highind == 1: # red is highest
if B >= G: # there is more blue than green
return 0 # max red and more blue than green shouldn't happen
# wavelength linearly changes from 645 to 580 as green increases
RGBlambda = 645 - RGBnorm[1] * (645. - 580.)

elif highind == 2: # green is max, range is 510 - 579
if R > B: # range is 510 - 579
RGBlambda = 510 + RGBnorm[0] * (580 - 510)
else: # blue is higher than red, range is 490 - 510
RGBlambda = 510 - RGBnorm[2] * (510 - 490)

elif highind == 3: # blue is max, range is 380 - 490
if G > R: # range is 440 - 490
RGBlambda = RGBnorm[1] * (490 - 440) + 440
else: # there is more red than green, range is 380 - 440
RGBlambda = 440 - RGBnorm[0] * (440 - 380)

return RGBlambda


And here is an accuracy check made on my older image of mercury lamp spectrum:
Accuracy check (synthetic test): x-axis is expected value, y-axis is algorithm output. Circles are the data, the ideal is presented by red line (y = x).

Please note that this was NOT produced by using a real photo of a spectrum. As someone who has some spectroscopic background I strongly discourage you from using this algorithm to get wavelength information. This was derived from a simplified version of simple wavelength-to-RGB algorithm. Nevertheless, I intend to make some photos of monochromatic light and test it ;-)

Sunday, February 21, 2010

How to make a .wav file with Python, revisited

I've noticed that my two-year old post about making a .wav file with Python is being read a lot so I forced myself to modify the code a little. Now it should be nicer to read and modify further according to your needs/wishes. I've also incorporated the comments of helpful readers Mike Axiak and Fingon. Thanks guys ;-)
The link to the source of this idea, code by Andrea Valle, does not work as I'm typing this so I can only thank him by saying his name. My contribution to all this is minor.

import numpy as N
import wave

def get_signal_data(frequency=440, duration=1, volume=32768, samplerate=44100):
"""Outputs a numpy array of intensities"""
samples = duration * samplerate
period = samplerate / float(frequency)
omega = N.pi * 2 / period
t = N.arange(samples, dtype=N.float)
y = volume * N.sin(t * omega)
return y

def numpy2string(y):
"""Expects a numpy vector of numbers, outputs a string"""
signal = "".join((wave.struct.pack('h', item) for item in y))
# this formats data for wave library, 'h' means data are formatted
# as short ints
return signal

class SoundFile:
def __init__(self, signal, filename, duration=1, samplerate=44100):
self.file = wave.open(filename, 'wb')
self.signal = signal
self.sr = samplerate
self.duration = duration

def write(self):
self.file.setparams((1, 2, self.sr, self.sr*self.duration, 'NONE', 'noncompressed'))
# setparams takes a tuple of:
# nchannels, sampwidth, framerate, nframes, comptype, compname
self.file.writeframes(self.signal)
self.file.close()

if __name__ == '__main__':
duration = 2
myfilename = 'test.wav'
data = get_signal_data(440, duration)
signal = numpy2string(data)
f = SoundFile(signal, myfilename, duration)
f.write()
print 'file written'