Determine the affine transform between two images or EBSD maps #mtexScript

imageAligner.m

function [M, U] = imageAligner(im0, varargin)
%%IMAGEALIGNER computes the affine transform between two images.
% [M, U]= IMAGEALIGNER(im0, im1) displays reference image im0 and deformed 
% image im1 in two seperate figures, and the user is asked to locate points
% appearing in both images (control points). Then, the transform matrix M 
% and the displacements U are given. Both images can be images as 2d array, 
% EBSD data or of class grain2d.
%
% It is assumed that coordinates in reference image X are transformed into 
% Y as follows:
%	Y = M*X + U
%
% [...]=IMAGEALIGNER(im0,'field0',...) plots the field 'field0' for
% reference image. im0 must be of class EBSD or grain2d.
%
% [...]=IMAGEALIGNER(...,im1,'field1') plots the field 'field1' for
% deformed image. im1 must be of class EBSD or grain2d.
%
% Example:
%	ebsd = EBSD.load('/Whatever/')
%	sem = imread('SEM.tif');
%	[M, U] = imageAligner(ebsd, 'orientations', sem);
%	ebsd_aligned=affinetrans(ebsd,M,U);		% Apply transform to EBSD data
%
%	% Superimpose aligned EBSD and image
%	plot(ebsd_aligned,ebsd_aligned.orientations,'micronbar','off')
%	hold on
%	im=imshow(sem);
%	im.AlphaData = 0.7;
%	 
% See also affinetrans
%

%% Dialog box
txt1='Click alternatively on both images to specify the locations of control points. At least 3 points are required.';
txt2='Right click or press any key when it is done.';
uiwait(msgbox({txt1,txt2},'Locations of control points','help'));

%% Process input arguments
if ~ischar(varargin{1})
	% imageAligner(data0, data1...
	field0=[];
	im1=varargin{1};
	if nargin==2
		%...)
		field1=[];
	else
		%...,field1)
		if ischar(varargin{2})
			field1=varargin{2};
		else
			error('The 3th argument must be the name of the field to be displayed')
		end		
	end
else
	% imageAligner(array, field0, data1...
	field0=varargin{1};
	im1=varargin{2};
	if nargin==3
		%...)
		field1=[];
	else
		%...,field1)
		if ischar(varargin{3})
			field1=varargin{3};
		else
			error('The 4th argument must be the name of the field to be displayed')
		end		
	end	
end

%% Open figure for transformed image
h1=figure('Name','Deformed data','NumberTitle','off');
if (isa(im1,'grain2d') || isa(im1,'EBSD'))
	if isempty(field1)
		plot(im1,'micronbar','off');
	else
		plot(im1, im1.(field1), 'micronbar', 'off');
	end
else
	imshow(im1)
end
hold on
s1=scatter([],[],[],'filled','MarkerEdgeColor',[0 0 0]);

%% Open figure for reference image
h0=figure('Name','Reference image','NumberTitle','off');
if (isa(im0,'grain2d') || isa(im0,'EBSD'))
	if isempty(field0)
		plot(im0,'micronbar','off');
	else
		plot(im0, im0.(field0), 'micronbar', 'off');
	end
else
	imshow(im0)
end
hold on
s0=scatter([],[],[],'filled','MarkerEdgeColor',[0 0 0]);

%% User specifies the control points
n=10;
Pts0=NaN(n,2);
Pts1=NaN(n,2);
i=1;
while true
	figure(h0);
	[x,y,b]=ginput(1);
	if isempty(x) || b~=1
		break
	end
	Pts0(i,:)=[x,y];
	set(s0,'XData',Pts0(1:i,1),'YData',Pts0(1:i,2),'CData',lines(i));
	figure(h1);
	[x,y,b]=ginput(1);
	if isempty(x) || b~=1
		break
	end
	Pts1(i,:)=[x,y];
	set(s1,'XData',Pts1(1:i,1),'YData',Pts1(1:i,2),'CData',lines(i));
	i=i+1;
end
n=i-1;
if n<3
	error('At least 3 control points must be defined.')
end
Pts0=Pts0(1:n,:);
Pts1=Pts1(1:n,:);

%% Compute M and U
% Centroids and local CSs
G0=mean(Pts0);
G1=mean(Pts1);
loc_coords_0=Pts0-repmat(G0,[n 1]);
loc_coords_1=Pts1-repmat(G1,[n 1]);

% Pseudo inverse
M=(loc_coords_0\loc_coords_1)';
U=G1'-M*G0';

end