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