% dip07_3.m	   gr24mar15    31oct17
% interp_sep  Linear Image Interpolation (coop. Philips Research)
%
clear all;
return

%%
iptsetpref('ImshowInitialMagnification',100); %%%%%%% note

a = double(imread('.\img\lena.tif'));	a=a(151:306,201:356);
% a = double(imread('.\img\cameraman.tif'));	a=a(51:206,80:181);

rep = 2;    L = 2.^rep;    % interp. factor (2 or 4)

b0 = interp2(a,rep,'nearest');	% n.neighbor interp. by 2, repeated rep times
b1 = interp2(a,rep,'linear');	% bilinear interp. by 2, repeated rep times
b3 = interp2(a,rep,'spline');	% cubic spline interp. by 2, repeated rep times
% figure(1); imshow(b0,[]);
% figure(2); imshow(b1,[]);
% figure(3); imshow(b3,[]);
figure(1); imshow([b0 b1 b3],[]);


%%
b3 = [b3  b3(:,end); [b3(end,:) b3(end,end)] ];  % add missing last row and column

N = 20;   % order of anti-imaging filter (N ~ 10*rep is needed if rep>1)
w = L * fir1(N, 1/L);	% lowpass, ome_c = pi/L (Hamming window) [try e.g. .7/L, 1/L, 1.3/L]
                        % aliasing is introduced if larger ome_c is used

ch = upsample(a',L)';    % cv = upsample(a,L);
CH = mean(abs(fft(ch,[],2))); % average 1-D horiz. spectrum of upsampled img
len = length(CH);		% note symmetry wrt pi/L
W = abs(fft(w,len));
	figure(1); ome=(0:len-1)/len*2;  % subplot(211); 
	semilogy(ome,CH/10^4,ome,W); xlabel('ome/pi'); drawnow;
%	subplot(212);	plot(cumsum(w));  % step response of filter

[nr,nc] = size(a);      % apply anti-imaging filter
chh = [];
for i = 1:nr,
   chh = [chh; conv(ch(i,:),w)];
	end;
c = [];
chhu = upsample(chh,L);
for i = 1:L*nc+N,
   c = [c conv(chhu(:,i),w)];
	end;
c = c(N/2+1:end-N/2,N/2+1:end-N/2);
figure(2);  imshow(c,[0 255]);

figure(3);  % compare bicubic interp with upsampling+filter
stem([upsample(a(128,1:50),L)]); hold on;  % available data
plot(c(128*L,1:50*L),'k');  %  upsample+filter
plot(b3(128*L,1:50*L),'r');  % cubic spline
hold off; axis([20,90,0,250]);

iptsetpref('ImshowInitialMagnification','fit');  %%%%%%%% note

return;

