% myisoper.m
function w=myisoper()
% USES: myplot.m, mygetpts.m
clear;
close all;
disp('Do not make x(n+1)=x(1)')
pts=mygetpts([-2 2 -2 2],1);
clf
x=pts(:,1)';
y=pts(:,2)';
n=length(x);    % Vertices of polygon
x(n+1)=x(1);    % Make the polygon closed
y(n+1)=y(1);
t(1) =0;        % Parameter values: length of polygon
for i=1:n,
  t(i+1)=t(i)+sqrt((x(i)-x(i+1))^2 + (y(i)-y(i+1))^2);
end;
for i=1:n,      % Refined division td
  for j=1:10,
    td(10*(i-1)+j) = t(i)*(9-j+1)/9.0 +t(i+1)*(j-1)/9.0;
  end;
end;
xstart= (x(2)-x(n))/sqrt((x(2)-x(n))^2 + (y(2)-y(n))^2);
xstop= xstart;
ystart= (y(2)-y(n))/sqrt((x(2)-x(n))^2 + (y(2)-y(n))^2);
ystop= ystart;
xd= spline(t,[xstart x xstop],td);  % Spline values at td
yd= spline(t,[ystart y ystop],td);
myplot(x,y,'k.',xd,yd,'MarkerSize',30), 
title('Spline curve through points');
axis equal, grid, pause;
m=length(xd);
%a1=polyare2([xd(1:m-1)' yd(1:m-1)']);  
%disp(['polyare2 = ' num2str(a1)])
a2=polyarea(xd,yd);
disp(['polyarea = ' num2str(a2)])
patch(xd,yd,'y')
xlabel(['Area of shaded region = ' num2str(a2) ' (MYSPL)'])
le=sum( ((xd(1:m-1)-xd(2:m)).^2+(yd(1:m-1)-yd(2:m)).^2).^(0.5) );
ylabel(['Length of boundary = ' num2str(le)])
     fprintf('sqrt(Area)/length= %12.4e\n',a2/(le*le) )
w= a2/(le*le);

function plt=myplot(varargin)
%  Purpose:  to label and plot in preferred font size etc. 
axes('FontSize',[14],'FontWeight','bold'); 
hold on;
plt0=plot(varargin{:});
set(plt0,'LineWidth',2);
hold off
% FILE  myplot.m  ends. 

function w=mygetpts(bds,option)
% MYGETPTS plots a grid and prompts the user to
%        enter by mouse click points (x,y) 
%        with bds(1) <x <bds(2), bds(3)<y<bds(4)
% nrx = number of grid subdivisions in x-direction
% nry = number of grid subdivisions in y-direction

if (nargin ==1)
  option=1;
end;

% First plot grid
nrx=10; nry=10;
figure
hold on;
axis([bds(1)-0.1 bds(2)+0.1 bds(3)-0.1 bds(4)+0.1]);
axis('equal')
for i=1:(nrx+1)
  x=(bds(2)- bds(1))*(i-1)/nrx+ bds(1);
  line([ x x],[ bds(3 )  bds(4)]);
end;
for j=1:(nry+1)
  y=(bds(4)- bds(3))*(j-1)/nry+ bds(3);
  line([bds(1) bds(2)],[ y y]);
end
n=0;
% Get points with mouse:
disp('Left mouse button picks points.')
disp('Right mouse button picks last point.')
but = 1;
while (but == 1)
  [xi,yi,but] = ginput(1);
  plot(xi,yi,'go')
  n = n+1;
  x(n,1) = xi;
  y(n,1) = yi;
end
% Plot the points
pts=[x y];
plot(pts(:,1),pts(:,2),'ko','LineWidth',1.5)
plot(pts(:,1),pts(:,2),'b-','LineWidth',1.5)
% Prompt for correction:
[d1,d2]=size(pts);
if (option)
  yn=input('Enter 1/0 to edit/not to edit points: ')
else
  yn=0;
end;
while ((yn==1)| (isempty(yn)==1))
  disp('Click with mouse the point to be edited. ')
  pointed=ginput(1);
  %find index of pointed point:
  index=1; tst=10;
  for i=1:d1
    tmp=norm(pointed-pts(i,:));
    if (tmp<tst)
      index=i; tst=tmp;
    end;
  end
  disp('Enter a new point: ');
  new=ginput(1);
  pts(index,:)=new;
% plot the grid again
  hold off;
  clf;
  axis([bds(1)-0.1 bds(2)+0.1 bds(3)-0.1 bds(4)+0.1]);
  axis('equal')
  for i=1:(nrx+1)
    x=(bds(2)- bds(1))*(i-1)/nrx+ bds(1);
    line([ x x],[ bds(3 )  bds(4)]);
  end;
  for j=1:(nry+1)
    y=(bds(4)- bds(3))*(j-1)/nry+ bds(3);
    line([bds(1) bds(2)],[ y y]);
  end
  hold on;
% plot the corrected points
  plot(pts(:,1),pts(:,2),'ko','LineWidth',1.5)
  plot(pts(:,1),pts(:,2),'b-','LineWidth',1.5)  
  if (option)
    yn=input('Enter 1/0 to edit/not to edit points: ');
  else
    yn=0;
  end;
end;
% End of correction loop
hold on;
axis([bds(1)-0.1 bds(2)+0.1 bds(3)-0.1 bds(4)+0.1]);
axis('equal')
for i=1:(nrx+1)
  x=(bds(2)- bds(1))*(i-1)/nrx+ bds(1);
  line([ x x],[ bds(3 )  bds(4)]);
end;
for j=1:(nry+1)
  y=(bds(4)- bds(3))*(j-1)/nry+ bds(3);
  line([bds(1) bds(2)],[ y y]);
end
plot(pts(:,1),pts(:,2),'ko','LineWidth',1.5)
plot(pts(:,1),pts(:,2),'b-','LineWidth',1.5)
% Save points
if (option)
  ok=input('Enter 1/0 to save/not save points: ');
else
  ok=1;
end;
if (ok==1)
  fid=fopen('mygetpts.dat','w');
  fprintf(fid, '%12.6f %12.6f \n',pts');
  fclose(fid);
  disp('Points saved in the file mygetpts.dat');
end;
w=pts;
% FILE mygetpts.m ends.

% FILE  myisoper.m  ends.
