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.

Link to Matlab script

Link to spectrum 1

Link to spectrum 2

Link to spectrum 3

Link to peak data

Link to Princeton Instruments data

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

  1. download all the files and pictures below into the same folder
  2. run Master.m
  3. when prompted, click on one ofthe 9 .bmp image files
  4. 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

Email

1) you can do anything 2) def. doing: when work rate X period > threshold 3) work rate = f(approach).

Dan Gareau Dan Gareau

Log In