0001 %WAVELENGTH_CALIBRATE - Takes a horizontal slice of a reference spectrograph, 0002 % and through user interaction finds the pixel locations 0003 % that correspond to known wavelengths which are then used 0004 % to calculate and output a calibrated wavelength scale. 0005 % 0006 % It's not necessary to use every row of a spectrograph to derive the 0007 % spectrogram. Although every column is always used. To specify a horizontal 0008 % slice (region), use the optional arguments y0 and h. 0009 % 0010 % This function uses the mercury peaks found in a compact fluorescent lamp's 0011 % (CFL) spectrum as a wavelength reference. It is possible to use other 0012 % reference peaks (e.g. neon) with this function by slightly modifying it. 0013 % 0014 % Syntax: lambda = wavelength_calibrate(img, roi) 0015 % 0016 % Inputs: 0017 % img - an image matrix or an image's filename. 0018 % roi - region of interest of spectrograph to process. 0019 % 0020 % Outputs: 0021 % lambda - a calibrated wavelength scale in nanometers. 0022 % 0023 % Example: 0024 % lambda = wavelength_calibrate('path/filename'); 0025 % h=100; y0=(size(I,1)-h)/2; lambda = wavelength_calibrate(I, [y0, h]); 0026 % 0027 % Other m-files required: image2spectrum 0028 % Subfunctions: none 0029 % MAT-files required: none 0030 % Other files required: CFL_Hg_peaks_labeled.png 0031 % 0032 % See also: IMAGE2SPECTRUM, POLYFIT 0033 % 0034 % Author: Jonathan Thomson 0035 % Work: 0036 % email: 0037 % Website: http://jethomson.wordpress.com 0038 % 0039 0040 function lambda = wavelength_calibrate(img, roi) 0041 0042 if (nargin < 1 || nargin > 2) 0043 usage('wavelength_calibrate(img, roi)'); 0044 end 0045 0046 % Look at CFL_Hg_peaks_labeled.png to see the numbered peaks and their 0047 % wavelengths. These are mercury (Hg) spectrum peaks. The zeros in 0048 % REF_PEAK_WAVELENGTHS are for non-mercury peaks. Note that your CFL's 0049 % spectrogram may show a second weaker peak near 407.783 nm next to peak 1 0050 % at 404.656 nm. This is a mercury peak, but is not shown in the image 0051 % CFL_Hg_peaks_labeled.png nor is it used for calibration. 0052 0053 %NIST Atomic Spectra Database 0054 %http://physics.nist.gov/PhysRefData/ASD/lines_form.html 0055 %REF_PEAK_WAVELENGTHS = [404.6565, 435.8335, 0, 0, 546.0750, 0, 579.0670]; 0056 0057 %NIST Handbook of Basic Atomic Spectroscopic Data 0058 %http://physics.nist.gov/PhysRefData/Handbook/Tables/mercurytable2.htm 0059 %"This handbook is designed to provide a selection of the most important 0060 %and frequently used atomic spectroscopic data in an easily accessible 0061 %format." 0062 REF_PEAK_WAVELENGTHS = [404.6563, 435.8328, 0, 0, 546.0735, 0, 579.0663]; 0063 0064 % This is a user defined constant, but peak numbers 2, 5, and 7 should 0065 % work fine if you are using a CFL as the reference. 0066 PEAKNUM = [2 5 7]; 0067 0068 scrsz = get(0, 'ScreenSize'); 0069 0070 if (~exist('roi','var') || isempty(roi)) 0071 roi = [0, 0]; 0072 end 0073 0074 Z = image2spectrum(img, 'rgb', roi); 0075 0076 refpknm = REF_PEAK_WAVELENGTHS(PEAKNUM); 0077 0078 %if rng is large and peak 7 is too close to peak 6, peak 6 might be 0079 %incorrectly picked. 0080 %rng = -10:10; 0081 rng = -5:5; 0082 0083 dataisgood = false; 0084 while (dataisgood == false) 0085 n = figure('Position', [1 1 0.9*scrsz(3) 0.9*scrsz(4)]); 0086 hold on 0087 plot(fliplr(Z)) 0088 m = max(max(Z)); 0089 a = axis; 0090 axis([a(1) a(2) 0 1.10*m]) 0091 0092 for li = 1:length(PEAKNUM) 0093 title(['Click on peak ' num2str(PEAKNUM(li)) '.']); 0094 [x, ignored] = ginput(1); % 2nd output of ginput() is mandatory 0095 x = round(x); 0096 0097 % The user cannot be expected to click exactly at the peak 0098 % location so we should search on either side of gpxl 0099 % (within a hardcoded range) to find the exact peak location. 0100 dx = x + rng; 0101 0102 [peak_rgb, pkloc] = max(Z(dx,:)); 0103 [peak, chan] = max(peak_rgb); % which channel has the highest peak 0104 0105 % x-coordinate/column number of the peak 0106 pkx(li) = dx(1) + pkloc(chan) - 1; 0107 0108 hold on, plot(pkx(li), peak, 'k*') 0109 end 0110 0111 title(['The reference peaks'' locations have been recorded. ' ... 0112 'Please answer the question in the console.']) 0113 rsp = input(['Did you pick the correct peaks? ' ... 0114 'If not, answer no to try again. [(y)es/(n)o/(q)uit]: '], 's'); 0115 if (~isempty(rsp) && lower(rsp(1)) == 'y') 0116 dataisgood = true; 0117 elseif (~isempty(rsp) && lower(rsp(1)) == 'q') 0118 close(n) 0119 error('wavelength_calibrate: instructed to quit by user.'); 0120 end 0121 0122 close(n) 0123 end 0124 0125 if (dataisgood == true) 0126 [p, s] = polyfit(pkx, refpknm, 2); 0127 lambda = polyval(p, (1:size(Z,1)).'); 0128 end 0129 0130 end