warning off

% phantom images are 128x128
s = 128; p = s*s;

% set determines which set of images is used in loadims
% shapes2 phantom
%     1: range of shifts 0.25:0.25:2
%     2: range of rotations 0.25:0.25:2
%     3: range of rotations 0.25:0.25:2 plus a constant shift of 2
% shepp-logan phantom
%     4: range of shifts 0.25:0.25:2
%     5: range of rotations 0.25:0.25:2
%     6: range of rotations 0.25:0.25:2 plus a constant shift of 2
% shapes2 phantom: large motion
%     7: range of shifts 1:12
%     8: range of rotations 1:12
% shepp-logan phantom: large motion
%     9: range of shifts 1:12
%     10: range of rotations 1:12

if (set<1 | set>10) error('need to declare set'); end;

% additive noise mean & stdev
loadims

rsoln = zeros(3,nim); rerr = zeros(1,nim); cr = zeros(1,nim);
asoln = zeros(3,nim); aerr = zeros(1,nim); ca = zeros(1,nim);

for i=1:nim,
	imt = Images(:,:,i);
	[cr(i) Mr] = rigid2Dreg(im0,imt);
	rsoln(1,i) = -Mr(1,3); rsoln(2,i) = -Mr(2,3); rsoln(3,i) = Mr(2,1);
	imr = imageWarp(imt,Mr);
	rerr(i) = sqrt(sum(sum( (imr-im0).^2 ))/p);

	[ca(i) Ma] = affine2Dreg(im0,imt);
	Mx = Ma(1,1); My = Ma(2,2); Sx = Ma(1,2); Sy = Ma(2,1);
	theta = abs(0.5*(- asin(Sx) + asin(Sy)));
	asoln(1,i) = -Ma(1,3); asoln(2,i) = -Ma(2,3); asoln(3,i) = theta;
	ima = imageWarp(imt,Ma);
	aerr(i) = sqrt(sum(sum( (ima-im0).^2 ))/p);

end;

rsoln(3,:) = rsoln(3,:)*180/pi; asoln(3,:) = asoln(3,:)*180/pi;

d = 1:nim;

fprintf('\nAverage condition number for rigid transform: %f\n',mean(cr));
fprintf('\nAverage condition number for affine transform: %g\n',mean(ca));

figure;
subplot(2,2,1); plot(d,rsoln(1,:),'b-',d,asoln(1,:),'r-.',d,actx,'k*');
mi = min([min(rsoln(1,:)),min(asoln(1,:)),min(actx)]); mi = mi - 0.1*abs(mi);
ma = max([max(rsoln(1,:)),max(asoln(1,:)),max(actx)]); ma = ma + 0.1*abs(ma);
if (mi==0), mi=-0.1; end; if (ma==0), ma = 0.1; end;
title('Translation along x'); axis([d(1)-1 d(nim)+1 mi ma])
ylabel('pixels');

subplot(2,2,2); plot(d,rsoln(2,:),'b-',d,asoln(2,:),'r-.',d,acty,'k*'); 
mi = min([min(rsoln(2,:)),min(asoln(2,:)),min(acty)]); mi = mi - 0.1*abs(mi);
ma = max([max(rsoln(2,:)),max(asoln(2,:)),max(acty)]); ma = ma + 0.1*abs(ma);
if (mi==0), mi=-0.1; end; if (ma==0), ma = 0.1; end;
title('Translation along y'); axis([d(1)-1 d(nim)+1 mi ma]);
ylabel('pixels');

subplot(2,2,3); plot(d,rsoln(3,:),'b-',d,asoln(3,:),'r-.',d,actt,'k*');
mi = min([min(rsoln(3,:)),min(asoln(3,:)),min(actt)]); mi = mi - 0.1*abs(mi);
ma = max([max(rsoln(3,:)),max(asoln(3,:)),max(actt)]); ma = ma + 0.1*abs(ma);
if (mi==0), mi=-0.1; end; if (ma==0), ma = 0.1; end;
title('Rotation angle'); axis([d(1)-1 d(nim)+1 mi ma]);
ylabel('degrees');

subplot(2,2,4); plot(1:nim,rerr,'b-',1:nim,aerr,'r-.'); title('RMS error');
axis([d(1)-1  d(nim)+1 0.07 0.11])





