function [bands, bm, times,fpkt] = pbands(f,Krow,Kcol,Kthin)
% [bands, bm, times] = pbands(f,Krow,Kcol,Kthin)
times = zeros(12,6);
times(1,:) = clock;
[fcx,fce,fch,fcr,frx,fre,frh,frr] = procim(f,Krow,Kcol);  % get 1D parabs
times(2,:) = clock;
fr1 = ridge1(fch);    % put ridge labels on
times(3,:) = clock;
frt = thinr(fr1,Kthin);    % thin out noise (fewer than Kthin adjacent ridges
times(4,:) = clock;
frl = label(frt);          % label 4-connected ridge components
times(5,:) = clock;
Kcoef = parcoe2(Kcol,Krow);    % get coef matrix for 2D parab fit
times(6,:) = clock;
[fa,fb,fc,fd,fe] = parabg2(f,Krow,Kcol,Kcoef,fr1);   % find 2D parabs
times(7,:) = clock;
[fpks,fi] = parpks(fa,fb,fc,fd,fe);      % pull out parab peak locations
times(8,:) = clock;
fpkl = label(fpks);     % label 4-connected 2D parab peaks
times(9,:) = clock;
fpkt = thinpks(fpkl,fc,fd);    % thin down to 1 representative of component
times(10,:) = clock;
fb1 = bands1(fpkt,frl,f,fa,fb,fc,fd,fe,fcx,fch,fcr,frx,frh,frr,fr1); 
times(11,:) = clock;
[bands,bm] = bands3(fb1,fpkt,6);  % produces equivalent bands and image mask
times(12,:) = clock;
