% CMSC/AMSC 460  Fall 2007
% Fitting the elephant
% Dianne O'Leary  

% Here is a program which you can use to experiment with
% different interpolations to the outline of an elephant.

% Load the elephant picture stored in elephant.tif

figure(1)
[I,map] = imread('elephant.tif');
imagesc(I)
colormap(gray)

% Pick the interpolation points for the upper profile.
% Note that the x-coordinates must be in ascending order
% in order for the piecewise linear and the cubic spline fits
% to be correct.  (I really should check the input for errors here.)

% Because of the way imagesc displays pictures, the y coordinates
% need to be negated to make our interpolated elephants right-side up.

disp('Pick interpolation points for the upper profile of the elephant.')
[xup,yup] = ginput;
nup = length(xup);
yup = - yup;

% Pick the interpolation points for the lower profile, but reuse the
% end points from the upper.

disp('Pick interpolation points for the lower profile.')
disp('The two endpoints from the upper will be reused.')
[xlo,ylo] = ginput;
xlo = [xup(1);xlo;xup(nup)];
ylo = [yup(1);-ylo;yup(nup)];
nlo = length(xlo);

disp('Interpolation points for the upper profile:')
[xup,yup]

disp('Interpolation points for the lower profile:')
[xlo,ylo]


% Compute and evaluate the interpolating functions:

xplot = linspace(min(xlo),max(xlo),100);
xloe = xplot;
xupe = xplot;

figure(2)
kplot = 1;

npts = nlo + nup - 2;

% First, the polynomial

   
  polyPup = polyfit(xup,yup,nup-1);
  yup_poly = polyval(polyPup,xupe);
  yup_plot = polyval(polyPup,xplot);

  polyPlo = polyfit(xlo,ylo,nlo-1);
  ylo_poly = polyval(polyPlo,xloe);
  ylo_plot = polyval(polyPlo,xplot);
  
  
  subplot(2,2,kplot)
  plot(xplot,yup_plot,xplot,ylo_plot)
  hold on
  plot(xlo,ylo,'*',xup,yup,'o')
  axis([0 400 max(-1000,min(ylo_plot)),min(500,max(yup_plot))])
  pp = sprintf('Polynomial elephant, %d points',npts);
  title(pp)
  kplot = kplot + 1;
  

% Next, the piecewise linear.
% I chose to use the Matlab function interp1, but
% using the book's functions would be fine.

ylo_plot = interp1(xlo,ylo,xplot);
yup_plot = interp1(xup,yup,xplot);

ylo_lin = interp1(xlo,ylo,xloe);
yup_lin = interp1(xup,yup,xupe);


subplot(2,2,kplot)
plot(xplot,yup_plot,xplot,ylo_plot)
hold on
plot(xlo,ylo,'*',xup,yup,'o')
pl = sprintf('Piecewise linear elephant, %d points',npts);
title(pl)
kplot = kplot + 1;

% Finally, the cubic spline

coef_lo = spline(xup,yup);
coef_up = spline(xlo,ylo);

ylo_plot = ppval(coef_lo,xplot);
yup_plot = ppval(coef_up,xplot);

ylo_spline = ppval(coef_lo,xloe);
yup_spline = ppval(coef_up,xupe);


subplot(2,2,kplot)
plot(xplot,yup_plot,xplot,ylo_plot)
hold on
plot(xlo,ylo,'*',xup,yup,'o')
cs = sprintf('Cubic spline elephant, %d points',npts);
title(cs)
kplot = kplot + 1;