% script to process silhouette
% first read scanned image of curve
a=imread('scan_6523211656_1.jpg');
b=a(230:530,158:558,1);
[I J]=find(b<150);
jmin=min(J);
jmax=max(J);
newI=zeros(1,jmax-jmin+1);
newJ=zeros(1,jmax-jmin+1);
count=1;
for ind=jmin:jmax
    vals=find(J==ind);
    if numel(vals) > 0
    newI(count)=round(sum(I(vals))/numel(vals));
    newJ(count)=ind;
    count=count+1;
    end
end

newJ=newJ(1:count-1);
newI=newI(1:count-1);
newI=300-newI;
plot(newJ,newI);
fid=fopen('silh.dat','w');
fprintf(fid,'%3d ',newJ);
fprintf(fid,'\n');
fprintf(fid,'%3d ',newI);
fclose(fid);