%parseCraterCountFromArcGISShapefile
%Edwin Kite
%1 March 2014

function [cx,cy,cr] = parseCraterCountFromArcGISShapefile(NameOfCraterCountArcGISShapefile)

a = shaperead(NameOfCraterCountArcGISShapefile);

picturesque = 0;


for cid = 1:max(size(a)) %crater id
    assert(sum(isnan(a(cid).X)+isnan(a(cid).Y))==2)
  
    xs = a(cid).X(isnan(a(cid).X + a(cid).Y)==0);
    ys = a(cid).Y(isnan(a(cid).X + a(cid).Y)==0);

    [xfit,yfit,Rfit] = circfit(xs,ys);
    if picturesque==1
        figure
        plot(xs,ys,'b.')
        hold on
        rectangle('position',[xfit-Rfit,yfit-Rfit,Rfit*2,Rfit*2],...
            'curvature',[1,1],'linestyle','-','edgecolor','r');
        title(sprintf('Best fit: R = %0.1f; Ctr = (%0.1f,%0.1f)',...
            Rfit,xfit,yfit));
        plot(xfit,yfit,'g.')
        xlim([xfit-Rfit-2,xfit+Rfit+2])
        ylim([yfit-Rfit-2,yfit+Rfit+2])
        axis equal
    end
    
    cx(cid) = xfit; cy(cid) = yfit; cr(cid) = Rfit;
end

