Matlab Code
Here is videotaped class I taught that covers the use of "f-mins" for curve fitting and a few other helpful tricks. It's not edited so it's long and boring but nice if these are the tricks you need! Check it out….
Spectrometer Calibration
This page contains all the necessary code to calibrate a Raman spectrometer. You will need a spectrometer capable of acquiring a spectrum that is a string of numbers, each the intensity at a particular wavelength and a calibration standard such as a Raman scattering target. The example given here uses a 785 nm LASER and Acetominophen standard but with a bit of modification, this technique should be able to work on emission standards such as Xenon or Argon lamps or work with any LASER wavelength.
Take a look at the Matlab code by clicking on this link.
First, download all the following files into one folder on your computer. The run the Matlab script "Run_n_show_result.m" This will show you the calibration on our system. Then you can replace the spectra with your own and modify the Matlab script accordingly to calibrate your spectrometer.
Abberations When Focusing Through an Optical Window
The programs contained in this folder simulate the aberrations caused focusing a beam through a beamsplitter. In most applications such as reflectance mode confocal microscopy, the beam is passed through the beam-splitter collimated such that there is no aberration. Sending focusing rays through however causes aberrations which are modeled here.
AbbCube.m simulates focusing through a 1-inch cube beam-splitter. The user is incouraged to first view the big picture by running the simulation as is but then set the zoom feature on line 27 to Ò1Ó to see what is happening at the focus where the aberration occurs.
AbbPlate.m simulates focusing through a thin pellicle beamsplitter at 45 degrees. The zoom settings are all at the very end of the program and the user is encouraged to play with them to see different scales in the simulation.
Source code in Matlab:
AbbCube.m
% AbbCube.m Dan Gareau 6-25-07
% This program finds the abberation caused by focusing through a
% thick optical window. The picture drawn will be the entire view
% unless the user sets the zoom feature (line 27) to 1.
clear all
close all
Nrays = 6;
f = 100 % mm focal length
d = 10; % mm lens diameter
NA = d/(2*f); % - numerical aperture of focus
n1 = 1; % - refractive index
n2 = 1.52; % - refractive index (window)
n3 = 1; % - refractive index
pos = 50; % mm front edge position of optical window
t = 25.4; % mm window thickness
shift = (n2-1)*t/n2; % P. 101 Ch. 4 prisms and mirrors
shift = t*(1-n1/n2); % mm shift of focus with window inserted
LA = t*sqrt((n3^2-NA^2)/(n2^2-NA^2))-n3/n2*t;
hh = d/2; % mm height of lens
zoom = 0; % zoom in to new focus to see the abberations
for i = 1:Nrays
h(i) = i/Nrays*hh;
figure(1)
plot([0 f],[h(i) 0],'k--')
hold on
plot([0 f],[-h(i) 0],'k--')
theta1(i) = atan(h(i)/f);
plot([0 pos],[h(i) (h(i) - tan(theta1(i))*pos)],'k-')
plot([0 pos],[-h(i) -(h(i) - tan(theta1(i))*pos)],'k-')
end
% plot the optical window front surface
plot([pos pos],[-hh hh],'k-')
plot([f f],[-hh/3 hh/3],'k--')
plot([(f+shift) (f+shift)],[-hh/3 hh/3],'k--')
plot([(f+shift-LA) (f+shift-LA)],[-hh/3 hh/3],'k--')
plot([f (f+shift)],[hh/4 hh/4],'k--')
text(f-(shift/4), hh/2.4,sprintf('s = %1.1f mm',shift))
% plot the optical window other surface
plot([(pos+t) (pos+t)],[-hh hh],'k-')
text(pos + t*0.2,hh*0.90,sprintf('n = %1.2f',n2),'fontsize',16)
text(pos - t*0.8,hh*0.90,sprintf('n = %1.2f',n1),'fontsize',16)
text(pos + t*1.2,hh*0.90,sprintf('n = %1.2f',n1),'fontsize',16)
% label the windothickness
plot([pos (pos+t)],0.7*[hh hh],'k--')
text(pos*1.05,0.75*hh,sprintf('t = %2.1f mm',t))
for i = 1:Nrays
Yinc(i) = h(i) - pos*tan(theta1(i)); % y position of beam incident on window
theta2(i) = asin(n1/n2*sin(theta1(i)));
Yexit(i) = Yinc(i) - t*tan(theta2(i)); % y position of beam exiting window
% plot the beam in the optical window n2
plot([pos (pos + t)],[Yinc(i) Yexit(i)],'k-')
plot([pos (pos + t)],-[Yinc(i) Yexit(i)],'k-')
c(i) = Yexit(i)/tan(theta1(i));
final(i) = pos + t + c(i);
% plot the beam exiting the wondow to focus
plot([(pos + t) final(i)],[Yexit(i) 0],'k-')
plot([(pos + t) final(i)],[-Yexit(i) 0],'k-')
end
xlabel('z [mm]','fontsize',18)
ylabel('y [mm]','fontsize',18)
if zoom == 1;
center = mean(final);
range = 1.3*(center - min(final));
axis([ (center - 2*range) (center + 2*range) -range*0.2 range*0.2])
plot([(f+shift-LA) (f+shift)],[range*0.16 range*0.16],'k--')
text(f+shift-LA*0.35, range*0.17, sprintf('LA = %1.0f µm', abs(LA*1000)))
end
set(gca,'fontsize',18)
AbbPlate.m
% Abberation.m
% dagonal window
clear all
close all
Nrays = 5;
f = 100; % mm focal length
d = 10; % mm lens diameter
NA = d/(2*f); % - Numerical Aperture
n1 = 1; % - refractive index
n2 = 1.51; % - refractive index (window)
n3 = 1; % - refractive index
pos = 50; % mm front edge position of optical window
t = 0.2; % mm window thickness
ang = 45; % deg angle of sheet beamsplitter
ang = pi*ang/180; % convert to radians
hh = d/2; % mm height of lens
dx = t*sin(ang);
dy = t*cos(ang);
V = (n2 - 1)/(n3 - n1); % Abbe V number P. 102 Ch 4 prisms & mirrors
ang0 = 0; % initial angle of window before slanting
% mm lateral displacement of axial beam
D = -t*(n2-1)/n2; % small angle approximation (not too good)
D = t*sin(ang0 - ang)/cos(ang); % mm lateral displacement of axial beam
D = t*sin(ang)*(1-cos(ang)/cos(ang0)); % mm lateral displacement of axial beam
% Draw Beamsplitter
figure(1)
plot([pos (pos + hh*tan(ang))],[0 hh],'k-')
hold on
plot([(pos - hh*tan(ang)) pos],[-hh 0],'k-')
Ht = t/sin(ang);
plot([(pos + Ht) (pos + hh*tan(ang) + Ht)],[0 hh],'k-')
plot([(pos - hh*tan(ang) + Ht) (pos + Ht)],[-hh 0],'k-')
xfinal = f*1.1;
for i = 1:(2*Nrays)+1
ii = i - Nrays;
if abs(ii) > 0
h(i) = ii*hh/Nrays;
theta1(i) = atan(h(i)/f);
y(i) = (h(i)/tan(theta1(i))-pos) / (1/tan(theta1(i))+ tan(ang));
x(i) = (h(i) - y(i))/tan(theta1(i));
plot([0 f],[h(i) 0],'k--')
plot([0 x(i)],[h(i) y(i)],'k-')
theta2(i) = pi/2 - ang - theta1(i); % angle of incidence to glass normal
theta3(i) = asin(n1/n2*sin(theta2(i))); % angle of transmission into window
l = t*tan(theta3(i));
ddx = t*tan(theta3(i))*cos(ang);
ddy = t*tan(theta3(i))*sin(ang);
y3(i) = y(i) - dy + ddy;
x3(i) = x(i) + dx + ddx;
plot([x(i) x3(i)],[y(i) y3(i)],'k-')
y4(i) = y3(i) - (xfinal-x3(i))*tan(theta1(i));
plot([x3(i) xfinal],[y3(i) y4(i)],'k-')
end
end
plot([0 xfinal],[D D],'k--')
zoom1 = 1; % zoom in to new focus to see the abberations
if zoom1 == 1;
deltafactor = (t*n2)*0.8;
deltafactory = (t*n2)*0.3;
ctrx = f ;
axis([ctrx-deltafactor ctrx+deltafactor -deltafactory 0.4*deltafactory ])
%plot([(f+shift-LA) (f+shift)],[range*0.16 range*0.16],'k--')
%text(f+shift-LA*0.35, range*0.17, sprintf('LA = %1.0f µm', LA*1000))
end
zoom2 = 0; % zoom in to new focus to see the abberations
if zoom2 == 1;
axis([100.0011 100.0014 -8e-4 -7.4e-4])
end
zoom3= 0; % zoom in to new focus to see the abberations
if zoom3 == 1;
axis([48 53 -4 4])
axis off
end
%axis([0 120 -6 6])
set(gca, 'fontsize',18)
xlabel('z [mm]')
ylabel('y [mm]')
%axis equal
The results from these simulations were used in the following:
Gareau D.S., Abeytunge S., Rajadhyaksha M., Full-pupil versus divided-pupil confocal line-scanners for reflectance imaging of human skin in vivo, Proc SPIE, 6443-31, 2007
Lens Design for Optical Window Focusing Correction
The program contained in this folder determines the correct launch angles for rays to focus through
an optical window. This might be a cover slip that holds a sample onto microscope slide. The code
works by creating a look up table for rays traced from the focus to the focusing optic, then sampling
the table to draw the forward rays. Master.m has many input variables to specify the geometry that
the user wishes to implement. The creation of the lookup table is done by the program called
CreateLookup.m and it is computationally demanding since many rays must be traced. The more rays
that are traced the better the program works in that the focus is never perfect but approximates focusing
to a point with more rays traced. Therefore, the user should set the variable Nrays to be as large as
possible and make sure that the focusing is acceptable. Using the current values, Nrays = 10,000 is
sufficient to keep the rays focusing to within 1 um of eachother and 100,000 rays gets the rays to
focus to within 0.1 um.
Happy focusing! --Dan
Source code in Matlab:
CreateLookup.m
% Abberation.m
hh = d/2; % mm height of lens
NA = n1*sin(atan(hh/f)) % [-] Numerical App.
y3 = zeros(Nrays,1);
th3 = zeros(Nrays,1);
outY = zeros(Nrays/10,1);
outTH = zeros(Nrays/10,1);
% trace from focus to lens surface
for i = 1:Nrays
th1 = asin(NA)*i/Nrays;% * 180/pi;
th2 = asin(n1/n2*sin(th1)); % Snell's baby!
th3(i) = asin(n2/n3*sin(th2)); % Snell's baby!
x1 = pos + t;
x2 = pos;
x3 = 0;
y1 = tan(th1) * (f-x1);
y2 = y1 + tan(th2)*t;
y3(i) = y2 + tan(th3(i))*pos;
if ifplot == 1
figure(1)
plot([x1 f],[y1 0],'r-')
hold on
plot([x1 x2],[y1 y2],'r-')
plot([x2 x3],[y2 y3(i)],'r-')
end
end
if ifplot == 1
xlabel('z [mm]','fontsize',18)
ylabel('y [mm]','fontsize',18)
set(gca,'fontsize',18)
end
Nper = Nrays/10;
for iii = 1:Nrays/Nper
i = iii*100;
outY(iii) = hh/(Nrays/Nper)*iii;
for ii = 1:Nrays
if y3(ii) < outY(iii)
outTH(iii) = th3(ii);
end
end
end
if ifplot == 1
figure(2)
plot(outY,outTH*180/pi,'ko')
xlabel('radius of launch from lens surface [mm]')
ylabel('angle of launch from lens surface [deg]')
end
Nrays = length(outTH);
save OUTlookup outY outTH f n1 n2 n3 pos t Nrays
PlotForward.m
n1f = n3; %Reverse refractive indices to go forward to focus
n2f = n2;
n3f = n1;
figure(3);clf
for i = 1:Nrays
th1 = outTH(i);
th2 = asin(n1f/n2f*sin(th1)); % Snell's baby!
th3 = asin(n2f/n3f*sin(th2)); % Snell's baby!
plot([0 pos],[outY(i) (outY(i) - tan(th1)*pos)],'k-')
hold on
plot([0 pos],[-outY(i) -(outY(i) - tan(th1)*pos)],'k-')
y1 = outY(i) - pos*tan(th1);
y2 = y1 - t*tan(th2);
y3 = 0;
x1 = pos;
x2 = pos + t;
x3 = pos + t + y2*tan(pi/2-th3);
plot([x1 x2],[y1 y2],'k--')
plot([x2 x3],[y2 y3],'k-')
plot([x1 x2],[-y1 -y2],'k--')
plot([x2 x3],[-y2 -y3],'k-')
end
plot([pos pos],[-max(outY) max(outY)],'k-')
plot([(pos+t) (pos+t)],[-max(outY) max(outY)],'k-')
xlabel('z [mm]','fontsize',18)
ylabel('y [mm]','fontsize',18)
set(gca,'fontsize',18)
Master.mk
% Master.m by Dan Gareau 6-25-07
% This program implements a ray trace using Snell's law to correctly focus
% light through an oprical window of high refractive index
clear all
close all
Nrays = 10000; % 10000 for 1 um artifical axial aberration, 100000 for 0.1um
f = 3; % mm focal length
d = 5.5; % mm lens diameter
n1 = 1.34; % - refractive index
n2 = 1.52; % - refractive index (window)
n3 = 1.34; % - refractive index
pos = 2; % mm front edge position of optical window
t = 0.5; % mm window thickness
ifplot = 1; % show the ray-trace from focus to obj lens
CreateLookup % yields OUTlookup.MAT --> outY outTH f n1 n2 n3 pos t Nrays
sprintf('... Done creating look up table for ray launch ...')
clear all
load OUTlookup
PlotForward
sprintf('... Done running the forward ray trace ...')
Feature Matching Mosaic Stitching
The program contained in this folder stitches a 3-by-3 mosaic together in Matlab
- download all the files and pictures below into the same folder
- run Master.m
- when prompted, click on one ofthe 9 .bmp image files
- when prompted, select the top left cluster (excluding one data point in the bottom right) do this by using the cross hairs to click once on the top left of the cluster ..... then click on the cluster's bottom right and hitting enter.
Source code in Matlab:
Master.m
% Master.m
% Dan Gareau 2008
% This Program builds a seemless mosaic from a batch of images by finding
% the relative positions via feature-matching.
if 1 == 1
clear all
close all
global frame1 frame2 InitialGuessXXX InitialGuessYYY MosaicSize Variation Overlap frameX frameY tempDir wind Showfit Xa Xb Ya Yb
%%%%%%%%%%%%%%%%%%%%%%%%% User Choices
InitialGuessXXX = 574; % [px] lateral displacement of sequential image pair X
InitialGuessYYY = 0; % [px] Y
wind = 15; % Size of the window around acceptable neighborhood
MosaicSize = 3; % The square mosic size is [MosaicSize] images wide
Showfit = 1; % 1 = graph the fitting algorithm in motion, 0 = don't (20% faster)
%%%%%%%%%%%%%%%%%%%%%%%%% Rock and Roll...
tic;
GetFrameFits % Fits for the x and y offesets for feature matching for all HORIZONTAL neighboring pair
toc
FlipEm % Flipps every even row of image offsets
disp(sprintf('Now running FlipEm.m which reverses the order of sequential rows'))
disp(sprintf('And Sniffer.m which detects statistically bad fits'))
OffXinit = OffsetX;
OffYinit = OffsetY;
showstats
sniffer % sniffs out the bad seeds, eronious fits
number = 0; % init
for i = 1:length(flags)
if flags(i)>0
number = number + 1; % Pretty straight forward
end
end
figure(1);clf
plot(OffXinit,OffYinit,'k*')
hold on
plot(OffsetX,OffsetY,'ro')
title('scatter plot of fits of image offsets')
legend('Original Fits','Corrected Fits')
xlabel('X offset [pix]')
ylabel('Y offset [pix]')
disp(sprintf('In this case, the number of eronious fits was %d out of %d',number,length(flags) ))
warning off
save output OffXinit OffYinit OffsetX OffsetY MosaicSize frameY frameX Showfit
clear all
close all
load output
tic;
makerows % calls the fitting results and FRAMES to assemble and save ROWS
toc
end % debug
clear all
close all
tic;
global showfit
Showfit = 1;
alignrows % Fits for the x and y offesets for feature matching for desending rows
toc;
clear all
close all
disp(sprintf('Now assembling the final mosaic'))
tic;
load output
load FinalResults
BuildMosaic % calls the fitting results and ROWS to assemble and save MOSAIC
toc
figure(2);clf
imagesc(mosaic)
colormap gray
GetFrameFits.m
% GetFrameFits.m
% Fits for the x and y offesets for featuremapping
% for all HORIZONTAL neighboring pairs
% Dan Gareau 2008
global InitialGuessXXX InitialGuessYYY frameX frameY tempDir Showfit
% Prompt User for Input Image Files
defImgLoadFolder = '';
[filename, tempDir] = uigetfile('*.bmp', ...
'Pick one of image files', fullfile(defImgLoadFolder,'*.bmp'));
figList = dir(fullfile(defImgLoadFolder, '*.bmp'));
% Sorting, first by the length of file name, then by alphanumerical order
len = [];
for i= 1 : 1 : length(figList)
len(i) = length( figList(i).name );
end
[len,idx] = sort(len);
for i = 1:length(figList)
figList2(i) = figList(idx(i));
end
sprintf('Now getting x-y coordinates for sequential frame fits')
k = 0;
ResultCounter = 0;
for j = 1:MosaicSize
kk = -1;
disp(sprintf('now fitting frame number %d of %d',(j-1)*MosaicSize+1,MosaicSize*MosaicSize))
for i = 1:MosaicSize
k = (j-1)*MosaicSize+i;
if mod(k,MosaicSize) > 0
ResultCounter = ResultCounter + 1;
frame1 = imread(fullfile(tempDir, figList2(k).name));
frame2 = imread(fullfile(tempDir, figList2(k+1).name));
if i == 1
[frameY frameX] = size(frame1);
end
if Showfit == 1
figure(2)
subplot(1,2,1)
imagesc(frame1)
axis equal
colormap gray
subplot(1,2,2)
imagesc(frame2)
axis equal
colormap gray
end
offx = -InitialGuessXXX; % initial guess at offset x
offy = InitialGuessYYY; % initial guess at offset y
if mod(j,2) == 0
offx = InitialGuessXXX; % initial guess at offset x
offy = InitialGuessYYY; % initial guess at offset y
end
start = [offy offx];
results = fminsearch('fit',start);
OffsetX(ResultCounter) = results(2);
OffsetY(ResultCounter) = results(1)*1e5; % VARIABLE FACTOR!!!
end
end
end
makerows.m
% makerows.m
% Dan Gareau 2008
global tempDir
ffigList = dir(fullfile(tempDir, '*.bmp'));
% Sorting, first by the length of file name, then by alphanumerical order
len = [];
for i= 1 : 1 : length(ffigList)
len(i) = length( ffigList(i).name );
end
[len,idx] = sort(len);
for i = 1:length(ffigList)
ffigList2(i) = ffigList(idx(i));
end
for rownum = 1:MosaicSize
clear prev deltaPOS cumsum Y1rel MaxFactorY newval SumFactorXX SumFactorX OffsetY OffsetX
clear OffsetNum index RowSize
load output
disp(sprintf('now building row number %d of %d',rownum,MosaicSize))
%index is the index of the beginning (leftmost) frame in the currrent
%row within the results file that has MosaicSize-1 frames per row
index(1,1) = (rownum-1)*(MosaicSize-1) + 1;
% SumfactorX is the difference between how wide the row is and how
% wide the row would be if the frames were just placed end to end
for iii = 1:MosaicSize-1
SumFactorXX(iii) = sum( OffsetX( index : index+iii-1 )- frameX );
end
SumFactorX = SumFactorXX(length(SumFactorXX));
% cumsum, the cumulative sum is the Y-offset between the LAST frame in
% the current row and the current frame
Icumsum = 0;
for iii = rownum*(MosaicSize-1):-1:index
Icumsum = Icumsum + 1;
cumsum(Icumsum,1) = -sum(OffsetY(iii:rownum*(MosaicSize-1)));
end
MaxFactorY = max(cumsum) - min(cumsum); % the biggest difference
if max(cumsum) < 0 % if last frame in row is LOWEST frame
MaxFactorY = min(cumsum);
end
if min(cumsum) > 0 % if last frame in row is HIGHEST frame
MaxFactorY = max(cumsum);
end
MaxFactorY = abs(MaxFactorY);
[A B] = max(cumsum); % A = y-offset max B = index of such
% deltaPOS = difference between top of LOWEST frame and top of 1st frame
C = length(cumsum) - B;
deltaPOS = 0;
for back = C:-1:1
deltaPOS = deltaPOS - OffsetY(back);
end
if cumsum(length(cumsum)) == max(cumsum) && max(cumsum) > 0
deltaPOS = 0;
end
if A < 0
deltaPOS = cumsum(length(cumsum));
end
Y1rel = MaxFactorY + deltaPOS; % first picture's y-position
clear merge1
merge1 = zeros(MaxFactorY+frameY+4, MosaicSize*frameX + SumFactorX+1) - 1;
[NNY NNX] = size(merge1);
for LocalFrameNum = 1:MosaicSize
if LocalFrameNum == 1
if cumsum(length(cumsum)) == min(cumsum) && cumsum(length(cumsum)) < 0
Y1rel = 0;
end
prev = Y1rel;
end
LocalOffsetNum = LocalFrameNum - 1;
OffsetNum = LocalFrameNum - 1 + (rownum-1)*MosaicSize-rownum+1;
if mod(rownum,2) == 0
k = LocalFrameNum + (rownum-1)*MosaicSize;
end
if mod(rownum,2) ~= 0
k = (rownum)*MosaicSize - LocalFrameNum + 1;
end
AddFrame = imread(fullfile(tempDir, ffigList2(k).name));
for j = 1:frameY
for i = 1:frameX
if LocalFrameNum == 1
ix = i;
iy = round(j + Y1rel);
else
ix = round(LocalOffsetNum*frameX + SumFactorXX(LocalOffsetNum)+ i);
iy = round( j + prev + OffsetY(OffsetNum));
end
if merge1(iy,ix) < 0
merge1(iy,ix) = AddFrame(j,i);
end
if merge1(iy,ix) >= 0
newval = mean([merge1(iy,ix) AddFrame(j,i)]);
merge1(iy,ix) = newval;
end
end
end
if LocalFrameNum > 1
prev = prev + OffsetY(OffsetNum);
end
end % end within row
if Showfit == 1
imagesc(merge1)
axis equal
colormap gray
drawnow
end
[tempY tempX] = size(merge1);
RowSize(rownum,1)=tempY;
RowSize(rownum,2)=tempX;
eval(['save row' num2str(rownum) ' merge1;']);
end % end of rows
save RowSize RowSize
alignrows.m
% alignrows.m
% Dan Gareau 12-2007
load output
global Showfit
for i = 1:MosaicSize
eval(['load row' num2str(i) ';']);
eval(['row' num2str(i) ' = merge1;']);
end
figure(1)
subplot(3,1,1)
imagesc(row1)
subplot(3,1,2)
imagesc(row2)
subplot(3,1,3)
imagesc(row3)
ResultCounter = 0;
for j = 1:MosaicSize-1
disp(sprintf('now fitting row number %d of %d',j+1,MosaicSize))
clear frame1
clear frame2
global frame1 frame2
ResultCounter = ResultCounter + 1;
eval(['frame1 = row' num2str(j) ';']);
eval(['frame2 = row' num2str(j+1) ';']);
[Ny Nx] = size(frame2);
if Showfit == 1
figure(2)
subplot(1,2,1)
imagesc(frame1)
axis equal
colormap gray
subplot(1,2,2)
imagesc(frame2)
axis equal
colormap gray
end
offx = -4; % -4 initial guess at offset x
offy = Ny*0.75;%0.75;%*0.7; % initial guess at offset y
start = [offy offx];
results = fminsearch('fitROWS',start);
OffsetXrows(ResultCounter) = results(2);% VARIABLE FACTOR!!!
OffsetYrows(ResultCounter) = results(1);
if j == MosaicSize-1
SizeLast = Ny;
end
end
save FinalResults OffsetXrows OffsetYrows SizeLast
BuildMosaic.m
% BuildMosaic.m
% Dan Gareau 2008
global Showfit
mosaic = zeros(frameY*MosaicSize,frameX*MosaicSize) - 1;
[Ny Nx] = size(mosaic);
for LocalFrameNum = 1:MosaicSize
eval(['load row' num2str(LocalFrameNum) ';']);
[Nfy Nfx] = size(merge1);
if LocalFrameNum == 1
iStart = Nx/2 - Nfx/2;
jStart = 0;
end
for j = 1:Nfy
for i = 1:Nfx
ix = round(iStart + i);
iy = round(jStart + j);
if mosaic(iy,ix) < 0
mosaic(iy,ix) = merge1(j,i);
end
if mosaic(iy,ix) >= 0 && merge1(j,i) > 0
newval = mean([mosaic(iy,ix) merge1(j,i)]);
mosaic(iy,ix) = newval;
end
end
end
if LocalFrameNum < MosaicSize
iStart = iStart + OffsetXrows(LocalFrameNum);
jStart = jStart + OffsetYrows(LocalFrameNum);
end
if Showfit == 1
figure(1);clf
imagesc(mosaic)
axis equal
colormap gray
drawnow
end
end % end within row
mosaic = uint8(mosaic);
imwrite(mosaic,'TaylorStitch.BMP','BMP')
fit.m
% fit.m
% Dan Gareau 2008
function err = fit(start)
global frame1 frame2 Showfit
[nny nnx] = size(frame1);
offy = round(start(1)*1e5); % VARIABLE FACTOR!!!
offx = round(start(2));
clear merge1
clear merge2
clear merged
merge1 = zeros(nny+abs(offy),nnx+abs(offx));
merged = merge1;
merge1 = merge1 - 1;
[nnny nnnx] = size(merge1);
merge2 = merge1;
if (offx > 0) && (offy > 0) % adding to the lower right
merge1(1:nny,1:nnx) = frame1;
merge2(offy+1:nnny,offx+1:nnnx) = frame2;
end
if (offx > 0) && (offy < 0) % adding to the upper right
merge1(abs(offy)+1:nnny,1:nnx) = frame1;
merge2(1:nny,offx+1:nnnx) = frame2;
end
if (offx < 0) && (offy > 0) % adding to the lower left
merge1(1:nny,abs(offx)+1:nnnx) = frame1;
merge2(offy+1:nnny,1:nnx) = frame2;
end
if (offx < 0) && (offy < 0) % adding to the upper left
merge1(abs(offy)+1:nnny,abs(offx)+1:nnnx) = frame1;
merge2(1:nny,1:nnx) = frame2;
end
if (offx > 0) && (offy == 0)
merge1(1:nnny,1:nnx) = frame1;
merge2(1:nny,offx+1:nnnx) = frame2;
end
if (offx < 0) && (offy == 0)
merge1(1:nny,abs(offx)+1:nnnx) = frame1;
merge2(1:nnny,1:nnx) = frame2;
end
err = 0; % initialize
count = 0;
for i = 1:nnnx
for j = 1:nnny
if (merge1(j,i) < 0) && (merge2(j,i) >= 0)
merged(j,i) = merge2(j,i);
end
if (merge1(j,i) >= 0) && (merge2(j,i) < 0)
merged(j,i) = merge1(j,i);
end
if (merge1(j,i) >= 0) && (merge2(j,i) >= 0)
merged(j,i) = (merge1(j,i)+merge2(j,i))/2;
err = err + abs(merge1(j,i)-merge2(j,i));
count = count + 1;
end
end
end
err = err/count;
if offx > nnx
err = err + 1e12;
end
if Showfit == 1
figure(3)
imagesc(merged)
axis equal
colormap gray
drawnow
end
fitROWS.m
% fitROWS.m
% Dan Gareau 2008
function err = fitROWS(start)
global frame1 frame2 Showfit
[nny1 nnx1] = size(frame1);
[nny2 nnx2] = size(frame2);
offy = round(start(1));
offx = round(start(2));%*1e5% VARIABLE FACTOR!!!
clear merge1
clear merge2
clear merged
if nnx1 >= nnx2
merge1 = zeros(nny2+offy,nnx1+abs(offx));
end
if nnx1 < nnx2
merge1 = zeros(nny2+offy,nnx2+abs(offx));
end
merged = merge1;
merge1 = merge1 - 1;
[nnny nnnx] = size(merge1);
merge2 = merge1;
if offx > 0 % adding to the lower right
merge1(1:nny1,1:nnx1) = frame1;
merge2(nnny-nny2+1:nnny,nnnx-nnx2+1:nnnx) = frame2;
end
if offx < 0 % adding to the lower left
merge1(1:nny1,nnnx-nnx1:nnnx-1) = frame1;
merge2(offy+1:nnny,1:nnx2) = frame2;
end
if offx == 0 % adding to the lower left
merge1(1:nny1,1:nnx1) = frame1;
merge2(offy+1:nnny,1:nnx2) = frame2;
end
err = 0; % initialize
count = 0;
for i = 1:nnnx
for j = 1:nnny
if (merge1(j,i) < 0) && (merge2(j,i) >= 0)
merged(j,i) = merge2(j,i);
end
if (merge1(j,i) >= 0) && (merge2(j,i) < 0)
merged(j,i) = merge1(j,i);
end
if (merge1(j,i) >= 0) && (merge2(j,i) >= 0)
merged(j,i) = (merge1(j,i)+merge2(j,i))/2;
err = err + abs(merge1(j,i)-merge2(j,i));
count = count + 1;
end
end
end
if count == 0
count = 1;
end
err = err/count;
if err == 0
err = 1e12;
end
if Showfit == 1
figure(3)
imagesc(merged)
axis equal
colormap gray
drawnow
end
FlipEm.m
% FlipEm.m
% Dan Gareau 2008
% Flipps every even row of image offsets
for i = 1:MosaicSize
if mod(i,2) ~= 0
iStart = 1+(MosaicSize-1)*(i-1);
iStop = iStart + MosaicSize - 2;
tempY = OffsetY(iStart:iStop);
tempX = OffsetX(iStart:iStop);
for ii = 1:MosaicSize-1
OffsetY(iStart+MosaicSize-1-ii) = -tempY(ii);
OffsetX(iStart+MosaicSize-1-ii) = -tempX(ii);
end
end
end
showstats.m
% showstats.m
global Xa Xb Ya Yb
figure(31);clf
plot(OffXinit,OffYinit,'k*')
hold on
%plot(OffsetX,OffsetY,'ro')
title('pick correct fit cluster region')
%legend('Original Fits','Corrected Fits')
xlabel('X offset [pix]')
ylabel('Y offset [pix]')
coord = ginput;
Xa = coord(1,1);
Xb = coord(2,1);
Ya = coord(2,2);
Yb = coord(1,2);
sniffer.m
% sniffer.m
% Dan Gareau 2008
global Xa Xb Ya Yb
Xa
Xb
Ya
Yb
global wind
MMMoffY = (Ya + Yb)/2;%median(OffsetY);
MMMoffX = (Xa + Xb)/2;%median(OffsetX);
STDoffY = (Yb-Ya)/2;%(std(OffsetY);
STDoffX = (Yb-Ya)/2;%std(OffsetX);
flags = zeros(length(OffsetY),1); % init, if outlier, flag = 1
% Detect outliers outside offset results neighborhood
for i = 1:length(OffsetY)
if OffsetY(i) > MMMoffY + wind | OffsetY(i) < MMMoffY - wind | OffsetX(i) > MMMoffX + wind | OffsetX(i) < MMMoffX - wind
%if OffsetY(i) > Yb | OffsetY(i) < Ya | OffsetX(i) > Xb | OffsetX(i) < Xa
flags(i) = 1;
end
end
resultcounter = 0;
summX = 0;
summY = 0;
for i = 1:length(OffsetY) % find mean of acceptable neighborhood
if flags(i) == 0
resultcounter = resultcounter + 1;
summX = summX + OffsetX(i);
summY = summY + OffsetY(i);
end
end
meanY = summY/resultcounter;
meanX = summX/resultcounter;
for i = 1:length(OffsetY)
if flags(i) == 1
OffsetY(i) = meanY;
OffsetX(i) = meanX;
end
end
