Skip to content

Instantly share code, notes, and snippets.

@Stoft1
Created May 18, 2015 17:23
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
Star You must be signed in to star a gist
Save Stoft1/2f3fd173f1ea3ccb2b45 to your computer and use it in GitHub Desktop.
Matlab example of CFL Calibration for PLab's SWB
% File: swbcfl.m
% Org: 20150218 Dave Stoft
% Revid: 20150323
% Descr: Example code to calculate pixel:wavelength mapping when pixel#
% for 436nm and 546nm peaks are known
%
% Concept:
% With an incidence angle to the grating of 45-deg, and the camera image
% plane thus parallel to the grating, the 524nm diffraction angle will be
% 0-deg; relative to the camera normal (perpendicular to the image plane
% and through the center of the lens). The 524nm pixel will be roughly at
% the mid-point 'X' of the camera image.
%
% As the wavelength extends below and above 524nm, the change in nm per
% pixel change will, itself, change -- it is not constant. However, the
% rate of change is linear.
%
% Note: When limited to the PLab spectrometer's useful range of 400-650nm,
% the correction per pixel is relatively small - on the order of 2nm. Is
% this still worthwhile to do? I think so, just because 1) it can be
% corrected, 2) that's what computers are good at and 3) it would seem to
% be totally transparent to the user.
%
% The spectrum will be spread either side of 524nm and the 436nm
% and 546nm CFL peaks will appear with 'X' pixel values below and above
% this 'center'. However, since the CFL has no peak at 524nm, the 436nm
% and 546nm peaks are used to calculate the 524nm pixel via simple trig
% functions. If the incedence angle is not exactly 45-deg, an error can
% result, but that error is much smaller than the amount of this
% correction so can probably be ignored.
%
% Given this configuration, the relationship between wavelength,
% diffraction angle and pixel 'X' position in the image plane can be
% calculated. Given only the pixel# of the CFL 436nm and 546nm peaks, a 1:1
% array can thus be generated which maps pixels to wavelength(nm).
%
% Assume:
% - Incidence angle (IA) from slit to diffraction grating is 45-deg
% - Camera normal is perpendicular to the diffraction grating plane
% - Visible spectrum is roughly centered within the camera image plane
% - The spectrum 'average' data has been taken and stored in an array
% - Peaks at 436nm (blue) and 546nm (upper green of doublet) are known
% - The pixel numbers of the max value of 436 and 546 peaks are known
%
% Calc:
% - With an IA of 45-deg, 524nm has a rel diffraction angle (DA) of 0-deg
% - For any pixel, we can calc the pixel 'x' offset to get the angle
% - Then, back-calc from the angle to get the wavelength
% - Using known 436 and 546 pixels, use trig to find pixel# of 524nm
% - Then, calc the wavlength for any pixel
% - So we get a 1:1 mapping of pixels to nm !
% ------------------------------------------------------------------------
function swbcfl( )
% ------------------------------------------------------------------------
% Declarations:
IA = 45.00; % From spectrometer construction
DGnm = 740; % Diffraction grating line spacing in nm
Zero_nm = 523.78; % Ideal center for 45.00 deg IA
Cal1_nm = 435.83; % Blue Hg line
Cal2_nm = 546.08; % Higher-nm peak of Grn Hg doublet
Pix436 = Pixel436; % Image pixel# of 436nm peak
Pix546 = Pixel546; % Image pixel# of 546nm peak
PixMax = 1600; % Total 'X' pix of the camera image
% ------------------------------------------------------------------------
% DEBUG: Read the REF image data from the file
[cfl_image, w_pix, h_pix] = read_image( [ 'P3U-CFL-03_TLapse.jpg'] );
% DEBUG: Read/average three line of colors from the image
y_line = 630;
for x_pix = 1:1600
rsum = 0;
gsum = 0;
bsum = 0;
for y_pix = -1:+1
rsum = rsum + int32(cfl_image(y_line+y_pix,x_pix,1));
gsum = gsum + int32(cfl_image(y_line+y_pix,x_pix,2));
bsum = bsum + int32(cfl_image(y_line+y_pix,x_pix,3));
end
pix(x_pix,1) = int16(double(rsum)/double(3));
pix(x_pix,2) = int16(double(gsum)/double(3));
pix(x_pix,3) = int16(double(bsum)/double(3));
cfl_uimage(x_pix) = ( pix(x_pix,1)+pix(x_pix,2)+pix(x_pix,3) ) / 3.0;
end
% ------------------------------------------------------------------------
% First, knowing the 436 and 546 pixels, find the 524nm pixel
%
% DVD
% | /|\
% | / | \
% | / | \
% FD / A | B \
% | / | \
% | / | \
% --------------- Image
% |DA436 | DA546|
% X436 X524 X546
%
% Notice that angles A and B share a common adjacent side
% Tan(A) = (X524-X436) / FD % A=abs(DA436), FD=focal distance
% Tan(B) = (X546-X524) / FD % B=abs(DA546)
% Since FD is in common
% (X524-X436)/Tan(A) = (X546-X524)/Tan(B)
% Then....
% Find the number of pixels between 436 and 546nm
% P = (Pix546-Pix436)
% So, X546 = P - X436
% And, by simple substitution...
% X436/Tan(A) = (P-X436)/Tan(B)
% Then, derive the pixel for 524nm...
% X436 = Tan(A)/Tan(B) * (P-X436)
% X436 = P*(Tan(A)/Tan(B)) - P436*(Tan(A)/Tan(B))
% X436 * (1+(Tan(A)/Tan(B)) = P*(Tan(A)/Tan(B))
% X436 = [ P*(Tan(A)/Tan(B)) ] / [ 1+(Tan(A)/Tan(B)) ]
% X546 = P - X436
% X524 = X436 + X436
% Next, derive the formula for diffraction angle
% Include the effect of 45-deg incidence angle
%
% Start with the basic diffraction equasion for non-zero IA
% d(sin(a) - sin(i)) = m*L
% where....
% d = grating line spacing
% a = grating exit angle
% i = incidence angle
% m = diffraction mode (1 here)
% L = lambda = wavelength of light
% Solving for angle (b)
% L / d = sin(a) - sin(i)
% sin(a) = (L/d) + sin(i)
% a = asin( (L/d) + sin(i) )
% Because of the polarity of the angle ...
% DA(nm) = asin( (Lnm/DGnm) - sin( IA*(2*pi/360) ) ) * (360/(2*pi));
% ========================================================================
% ========================================================================
% Now, do the calculations for real
% Diffraction angle for 524nm -- we've declared it to be zero
DA524 = 0.0;
% Diffraction angle for 436nm (trig functions work in radians)
DA436 = abs( asin( (Cal1_nm/DGnm) - sin(IA*(2*pi/360) ) ) * (360/(2*pi)) );
% Diffraction angle for 546nm
DA546 = abs( asin( (Cal2_nm/DGnm) - sin(IA*(2*pi/360) ) ) * (360/(2*pi)) );
% Find the image pixel span from Cal points: 436nm to 546nm
PixSpan = Pix546 - Pix436;
% Find the tangent ratio -- because the ZeroDA distance is common
TanAB = (tan(DA436*(2*pi/360))) / (tan(DA546*(2*pi/360)));
% Calc the #pix on the 436nm side -- the pix distance at image from 524nm
NumPix436 = int16( (PixSpan*TanAB) / (1+TanAB) );
% Now calc the ZeroDA distance which is in common
ZPixFD = double(NumPix436) / ( tan(DA436*(2*pi/360)) );
% Find the #pix from ZeroDA to 546nm image pix
NumPix546 = PixSpan - Pix436;
% Find the image pix for theoretical ZeroDA 524nm
Pix524 = Pix436 + NumPix436;
% Now, we have all we need to calc all pixle wavelengths
%
% Calc all Negative diffraction angles (per pixel)
f400 = 0; % debug use
for pix = Pix524:-1:1
% We now know focal distance 'ZPixFD' so we can find the diffraction
% angle for any pixel distance relative to 524nm
% Remember, for nm<ZeroDA, the diffraction angles are all negative
%
% First, find the DA angle relative to the ZeroDA (camera normal)
% L = atan( pixel_offset_from_524nm / FD )
DA(pix) = atan( double(Pix524-pix) / ZPixFD ) * (360/(2*pi));
% Then find the nm for this angle
% Start with the basic diffraction equasion for non-zero IA
% d(sin(a) - sin(i)) = m*L
% L = d(sin(a) - sine(i))
% NM(pix) = DGnm * ( sin(DA(pix)*(2*pi/360)) - sin( IA*(2*pi/360) )
PIXnm(pix) = DGnm * ...
( sin(IA*(2*pi/360)) - sin(DA(pix)*(2*pi/360)) );
% DEBUG: Identify the 400nm pixel
if ( f400 == 0 )
if ( PIXnm(pix) < 400 )
Pix400 = pix;
f400 = 1;
end
end
end
% Calc all Positive diffraction angles (per pixel)
f650 = 0;
for pix = (Pix524+1):PixMax
DA(pix) = atan( double(pix-Pix524) / ZPixFD ) * (360/(2*pi));
PIXnm(pix) = DGnm * ...
( sin(IA*(2*pi/360)) + sin(DA(pix)*(2*pi/360)) );
% DEBUG: Identify the 650nm pixel
if ( f650 == 0 )
if ( PIXnm(pix) > 650 )
Pix650 = pix;
f650 = 1;
end
end
end
% PIXnm(1:1600) now contains the (nm) values for every 'X' pixel in the
% spectrum. This is the end of the necessary calculations.
% ========================================================================
% ========================================================================
% ------------------------------------------------------------------------
% DEBUG: Find the near-linear nm/pix rate of change (look near 524nm)
Zslope = (PIXnm(Pix524+10)-PIXnm(Pix524-10)) / 20.0;
% DEBUG: Calc the offset from a linear estimate due to disperaion
for pix = Pix524:-1:1
% Linear estimate (nm) = 524nm - (nm/pix * pix_offset_from_524nm)
linear_est(pix) = double(PIXnm(Pix524)) - (Zslope * double(Pix524-pix));
% Find nm difference from linear est
nmdif(pix) = double(PIXnm(pix)) - linear_est(pix);
end
for pix = Pix524:1600
linear_est(pix) = double(PIXnm(Pix524)) + (Zslope * double(pix-Pix524));
nmdif(pix) = linear_est(pix) - double(PIXnm(pix));
end
pmin = Pix400;
pmax = Pix650;
% ------------------------------------------------------------------------
% DEBUG: Plot CFL and Dispersion on same plot
figure( 'position', [ 200 400 800 350] );
[AX,H1,H2] = ...
plotyy( (pmin:pmax), cfl_uimage(pmin:pmax), (pmin:pmax), nmdif(pmin:pmax) );
xlabel( 'Pixel Position', 'FontSize',12 );
set(get(AX(1),'Ylabel'),'String','Intensity (unitless)');
set(get(AX(2),'Ylabel'),'String','Dispersion(nm)');
set(get(AX(1),'Ylabel'),'FontSize',14);
set(get(AX(2),'Ylabel'),'FontSize',14);
%ylabel( 'Intensity (unitless)', 'FontSize',12 );
title( 'PLab3U CFL Spectra + Dispersion', 'FontSize', 12 );
set(AX(1),'xlim',[pmin pmax]);
set(AX(1),'ylim',[0 150]);
set(AX(2),'xlim',[pmin pmax]);
set(AX(2),'ylim',[-1 +2]);
grid on;
return;
% ------------------------------------------------------------------------
% DEBUG: Read the image data
function [new_image, w_pix, h_pix] = read_image( image_name )
% Fetch the image from the file
new_image = imread( image_name );
% Find the image dimensions
[ h_pix, w_pix, colors ] = size(new_image);
return;
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment