%try
 
  %uncomment for standalone operation
  fn='hori2Dgran_rough10_run1_050411_222455_0003'
%  fn='hori2Dgran_rough1_run21_040819_204811_0046'
%  fn='hori2Dgran_rough10_run50_050413_035951_0203';
%  fn='hori2Dgran_rough10_run56_050413_215724_0000';
  Nfac=1;
  %uncomment for standalone operation
  
%% Define Constants
  rev='006';
  diary([fn '.diary']);  % Turn on diary

  % Define default code behavior
  remove=1;   % remove bin file 1=yes 0=no
  plotit=1;   % plots during calculation 1=yes 0=no
  nplot=0;    % number of frames to skip before ploting
  pauseit=0;  % pause after each frame 1=yes 0=no
  saveit=1;   % save results 1=yes 0=no
  savenn=1;   % save nn file (Cummulitive histogram for each pixel) 1=yes 0=no
  refine=1;   % refine position 1=yes 0=no
  recalcbg=0; % force bg recalc 1=yes 0=no
  recalcDw=0; % force Dw recalc 1=yes 0=no
  recalc=0;   % force track recalc 1=yes 0=no
  Nstep=3;    % Number of processing steps
  
  varstr={'bg1','bg2','bgmx','bgmn','nmx','nmn','Dw','D','DD','w','ww','px','py','Npf','done','chi'};  % track variables to save
  varltr=varstr;                                                                      % track variables to load
  varsbg={'bg1','bg2','bgmx','bgmn','static','mod1','mod2','wmod','mn','mx'};         % bg variables to save
  varlbg={'bg1','bg2','bgmx','bgmn','static','wmod'};                                 % bg variables to load
  
  overide('track',rev,'_modcb',fn);  % override default Code Behavior
  
  % Default tracking parameters
  D=25;                    % first guess for D
  w=.8;                    % first guess for w
  MaxBalls=400;            % Maximum number of balls to find
  GoodBallsCutOff=5;       % Cutoff peak height for ball (not D independent)
  FracPart=.02;            % fraction of particle near edge to consider
  MinSepFrac=.65;          % Minimum separation between balls as a fraction of D
  clipx=1;                 % number of rows to clip at top
  clipy=1;                 % number of rows to clip at right
  mxcout=1.15;             % max value for scaled image
  mncout=-.15;             % min value for scaled image
  testfr=1;                % Frame(s) to use to get initial values.  Can be list.
  MaxPass=10;              % Maximum number of passes.

  overide('track',rev,'_modtp',fn);  % override default Tracking Parameter
  
  % Default background parameters
  
  graylevels=256;   % number of gray levels
  maskwid=30;       % width of mask to lookfor second peak
  mnmxfrac=.01;     % Nbg*mnmxfrac lowest and Nbg*(1-mnmxfrac) highest used to calculate mn and mx
  staticfrac=.02;   % Nbg*(1-staticfrac) must be static for pixel to be static

  overide('track',rev,'_modbg',fn);  % override default BackGround parameter

%% Check files and Setup
  fnimg=sprintf('%s.img',fn); % Image filename
  if(~exist(fnimg))           
    dalraw2img(fn,remove);    % Convert File if needed
  end

  [Nf info]=dalinfo(fnimg);     % get file info
  
  % image size
  Nx=info.Nx;            % Number of pixels in x-direction
  Ny=info.Ny;            % Number of pixels in y-direction

  Nbg=Nf;                % Number of frames to process for the background correction
  Nf=Nf;                 % Process only Nfac fraction of the frames
  frames=1:Nf*Nfac;      % List of Frames to process
  sfr=1;                 % current frame number
  efr=Nf*Nfac;           % last frame


  px=zeros(MaxBalls,Nf);   % x component of particle positions
  py=zeros(MaxBalls,Nf);   % y component of particle positions
  DD=zeros(1,Nf);          % Diameter for each frame
  ww=zeros(1,Nf);          % width for each frame
  Npf=zeros(1,Nf);         % Number of particles per frame
  chi=zeros(MaxPass,Nf);   % Goodness of fit for each frame
  done=zeros(1,Nf);        % List of progress
  step=0;                  % current step of process 0=none 1=bg, 2=track, 3=refine
  
  trfn=sprintf('%s%s_%s.mat','track',rev,fn);  % Track file
  bgfn=sprintf('%s%s_%s.mat','bg',rev,fn);     % background file
  Dwfn=sprintf('%s%s_%s.mat','Dw',rev,fn);     % Diameter file
  
  if (exist(trfn))
    load(trfn,'done');
    step=max(done);
    [jj sfr]=max(done<step); % Start from frame frames(sfr);
    if(sfr~=1)
      step=jj;
    end
    if (step>=2)
      load(trfn);
    end
  end

  fprintf(1,'processing file %s from step (%d/%d).\n',fnimg,step+1,Nstep)

  overide('track',rev,'_mod',fn);  % override default behavior
  
%% Calculate first try background (Step=1)
  if (step<1)
    % check to see if full background used
    if (exist(bgfn))
      load(bgfn,'wmod');
      if (max(wmod(:))<Nf/2)
        recalcbg=1;
        recalcDw=1;
      end
    end

    if (exist(bgfn) & recalcbg==0)  % Load background if previously calculated
      fprintf(1,'Loading Background Correction...');
      load(bgfn,varlbg{:});
      fprintf(1,'Done.\n');
    else
      if(recalcbg==0 & UpdateOldFiles('%s%03d_%s.mat','bg',[4 5],fn,bgfn,1)) % Lookfor older versions
        fprintf(1,'Done.\n');
        fprintf(1,'Loading Converted Background Correction...');
        load(bgfn,varlbg{:});
        fprintf(1,'Done.\n');
      else
        fprintf(1,'Computing Background Correction...');
        if exist(['nn_' bgfn])
          load(['nn_' bgfn],'nn','bb');
        else
          [nn bb]=himage4(fnimg,Nbg);   % Cummulitive histogram for each pixel
        end
        save(['nn_' bgfn],'nn','bb');
        [bgmx,bgmn,static,mod1,mod2,wmod,mn,mx]=...   % Calculate most probable max and min background
          bgcalc4(nn,bb,Nbg,Nx,Nx,graylevels,maskwid,mnmxfrac,staticfrac,mncout,mxcout);
        fprintf(1,'Done.\n');
        fprintf(1,'Saving Background Correction...');
        bg1=bgmx';
        bg2=bgmx'-bgmn';
        static=static';
        save(bgfn,varsbg{:});
        fprintf(1,'Done.\n');
        clear('nn');
        if (~savenn)
          fprintf(1,'deleting nn file...');
          delete(['nn_' bgfn]);
        end
        fprintf(1,'Done.\n');
      end
    end
    done=ones(1,Nf);        % List of progress (completed 1st step)
    step=1;
  end

%% Begin tracking

  if (step>=1)
    if (exist(trfn) & recalc==0)  % Load previous track file
      fprintf(1,'Loading track file...');
      load(trfn,varltr{:});
      fprintf(1,'Done.\n');
    elseif(recalc==0 & UpdateOldFiles('%s%03d_%s.mat','track',[4 5],fn,trfn,1)) % Lookfor older versions
      fprintf(1,'Done.\n');
      fprintf(1,'Loading Converted track file...');
      load(trfn,'px','py','Npf');
      fprintf(1,'Done.\n');
      nmx=0;
      nmn=0;
      if(UpdateOldFiles('%s%03d_%s.mat','Dw',[4 5],fn,Dwfn,1))
        load(Dwfn)
        fprintf(1,'Done.\n');
      end
      fprintf(1,'Saving Converted track file...');
      save(trfn,varstr{:});
      fprintf(1,'Done.\n');
    end

    if((step>=1) & (efr>sfr)) % First Pass to refine bg, Dw, and get initial px, py, Npf, chi

      [smip,ww3,mk,NNx,NNy,dD]=calcsmip(D,w,Nx-clipx,Ny-clipy,FracPart);  % calculate ideal particle based on current info
      ox=fix((dD-clipx)/2);                                               % position offset relative to raw image
      oy=fix((dD-clipy)/2);
      [xxbox yybox]=ndgrid(1:Nx-clipx,1:Ny-clipy);
      edgmask=zerofill(boxmask([(1+D/4)*[1 1] (Nx-clipx-D/4)*[1 1] 1+D/4 Ny-clipy-D/4],xxbox,yybox),NNx-Nx+clipx,NNy-Ny+clipy);
      rect=[oy+1 Ny+oy ox+1 Nx+ox];

      
      if(plotit==1)  % plot progress
        hs=simage(zeros(Nx,3*Ny));                
        hold on
        hp=plot(1:2,1:2,'k.');
        hold off;
      end
      
      if (max(done)==1)  % zero if first time through
        nmx=0;
        nmn=0;
        bgmn=0;
        bgmx=0;
      else
        bg1=bgmx./(nmx+eps);
        bg2=bg1-bgmn./(nmn+eps);
      end
      inx=(1:(Nx*Ny))';
      t0=clock;
      for nf=sfr:efr;
        raw=readimgarr(fnimg,Nx,Ny,frames(nf));    % load raw image
        craw=clip((bg1-raw)./bg2,mncout,mxcout);    % background correction and clip large values
        img=craw(clipx+1:end,clipy+1:end);          % clip crud at top and side;
        
        if (Npf(nf)==0) % Crude particle locator
          chiimg=1+(-2*conv2l(img,smip.^2,'full')+conv2l(img.^2,smip,'full'))./ww3;    % best fit ignoring overlap  
          [Npf(nf),spx,spy]=findpeaks(1./chiimg,mk,GoodBallsCutOff,MinSepFrac*D);      % find all maxima > CutOff with separation > Frac*D
          px(1:Npf(nf),nf)=spx;                                                        % store in px and py
          py(1:Npf(nf),nf)=spy;
        else            % or use previous location
          spx=px(1:Npf(nf),nf);  % store location in temp variable
          spy=py(1:Npf(nf),nf);
        end

        
        [cxy over]=pgrid(spx-ox,spy-oy,Nx-clipx,Ny-clipy,[1 Nx-clipx 1 Ny-clipy],Npf(nf),dD+2,0);  % create local grid centered on each particle and overlap matrix
        ci=ipf(abs(cxy),D,w);                                                % Calculate ideal image. Use ci=(1-tanh((abs(cr)-D/2)*w))/2; for small speedup.
        ri=img;
        di=(ci-ri);
        chi(done(nf),nf)=sum(di(:).^2);
        
        [dpx dpy]=cidp2(cxy,over,di,Npf(nf),D,w);
        spx=spx+dpx;
        spy=spy+dpy;
        [cxy over]=pgrid(spx-ox,spy-oy,Nx-clipx,Ny-clipy,[1 Nx-clipx 1 Ny-clipy],Npf(nf),dD+2,0);  % create local grid centered on each particle and overlap matrix
        ci=ipf(abs(cxy),D,w);                                                   % Calculate ideal image. Use ci=(1-tanh((abs(cr)-D/2)*w))/2; for small speedup.
        di2=(ci-ri);
        
        nr=0;
        delD=1;
        while((abs(delD)>.0001) & nr<10)
          [delD delw]=cidDw(abs(cxy),di2,D,w); D=D+delD;w=w+delw; ci=ipf(abs(cxy),D,w);di2=(ci-ri);
          fprintf('.');
          nr=nr+1;
        end
        Dw(nf,:)=[D w delD delw];
        
        [dpx dpy]=cidp2(cxy,over,di2,Npf(nf),D,w);
        spx=spx+dpx;
        spy=spy+dpy;
        [cxy over]=pgrid(spx-ox,spy-oy,Nx-clipx,Ny-clipy,[1 Nx-clipx 1 Ny-clipy],Npf(nf),dD+2,0);  % create local grid centered on each particle and overlap matrix
        ci=ipf(abs(cxy),D,w);                                                   % Calculate ideal image. Use ci=(1-tanh((abs(cr)-D/2)*w))/2; for small speedup.
        di2=(ci-ri);

        nr=0;
        delD=1;
        while((abs(delD)>.0001) & nr<10)
          [delD delw]=cidDw(abs(cxy),di2,D,w); D=D+delD;w=w+delw; ci=ipf(abs(cxy),D,w);di2=(ci-ri);
          fprintf('.');
          nr=nr+1;
        end
        Dw(nf,:)=[D w delD delw];
  
        [dpx dpy]=cidp2(cxy,over,di2,Npf(nf),D,w);
        spx=spx+dpx;
        spy=spy+dpy;
        [cr over]=pgrid(spx-ox,spy-oy,Nx,Ny,[1 Nx 1 Ny],Npf(nf),dD+2,1);  % create local grid centered on each particle and overlap matrix
        ci=ipf(cr(1:Nx-clipx,1:Ny-clipy),D,w);                                                   % Calculate ideal image. Use ci=(1-tanh((abs(cr)-D/2)*w))/2; for small speedup.
        di2=(ci-ri);

        mkmn=(cr<=D*(1-.1)/2);     
        mkmx=(cr>=D*(1+.1)/2);
        nmx=nmx+mkmx;
        nmn=nmn+mkmn;
        bgmx=bgmx+raw.*mkmx;
        bgmn=bgmn+raw.*mkmn;
        
        px(1:Npf(nf),nf)=spx;                                                        % store in px and py
        py(1:Npf(nf),nf)=spy;
        
        done(nf)=done(nf)+1;
        chi(done(nf),nf)=sum(di2(:).^2);

        if (plotit==1 & rem(nf-1,nplot+1)==0)
%          set(hs,'Cdata',[bgmx./(nmx+eps) raw bgmn./(nmn+eps)]);
          set(hs,'Cdata',[di ri di2]);
          caxis([-.4 .4]);
          set(hp,'xdata',spy-oy+Nx,'Ydata',spx-ox);
          drawnow;
        end
        fprintf(1,'Pass 1:: D=%6.3f w=%5.3f delD=%5.2e delw=%5.2e [%6.2f/%6.2f] fr=%04d et=%1.2f est.=%1.2f Np=%d id=%03d\n',Dw(nf,:),chi(done(nf)+(-1:0),nf),nf,etime(clock,t0)/60,etime(clock,t0)*((efr-sfr+1)/((nf-sfr+1))-1)/60,Npf(nf),str2num(info.id));
      end
      [step sfr]=max(done); 
      fprintf(1,'saving progress...step (%d/%d)...',step+1,Nstep)
      tpx=px;
      tpy=py;
      NBm=max(Npf(:));
      px(NBm+1:end,:)=0;
      py(NBm+1:end,:)=0;
      save(trfn,varstr{:});
      px=tpx;
      py=tpy;
      fprintf('Done\n');
    end
  end
  
%   
%   
%   
%   
%   
%   
%   
%   
%   
%   
%   
%   
%   
%   
%   
%   
%   
%   [dD,xxs,yys,smip,ww3,mk,NNx,NNy,sx,sy,rect,hs,hp,c1,wwm1,c2,wwm2,c3,edgmask]=DwdependentVar(D,w,Nx,Ny,static,refine);
%   [c1x c1y]=size(c1);
%   [c2x c2y]=size(c2);
% 
% 
%   Dwfn=sprintf('%s%s_%s.mat','Dw',rev,fn);
%   if (exist(Dwfn) & recalcDw==0)
%     fprintf(1,'Loading diameter and width...\n');
%     load(Dwfn,'D','w');
%     [dD,xxs,yys,smip,ww3,mk,NNx,NNy,sx,sy,rect,hs,hp,c1,wwm1,c2,wwm2,c3,edgmask]=DwdependentVar(D,w,Nx,Ny,static,refine);
%   else
%     if(recalcDw==0 & UpdateOldFiles('%s%03d_%s.mat','Dw',[4],fn,Dwfn,1))
%       fprintf(1,'Loading converted diameter and width...\n');
%       load(Dwfn,'D','w');
%     else
%       %calc position for 1st file
%       fprintf(1,'Caclulation diameter and width...\n');
%       %[junk nf]=max(find(polyok))
%       nf=1;
%       raw=readimgarr(fnimg,Nx+clipx,Ny+clipy,frames(nf));
%       raw=clip((bg1-raw)./bg2.*(1-static),mncout,mxcout);
%       img=raw(clipx+1:end,clipy+1:end);% clip crud at side
% 
%       ws=(1-(conv2l((conv2l(img>ithresh,c1,'full').*wwm1)>=1,c1,'same')<1)); % wall mask
%       ws=ws((c1x-1)/2+1:end-(c1x-1)/2,(c1y-1)/2+1:end-(c1y-1)/2);
%       wp=(1-(conv2l((conv2l((img>ithresh)-ws,c2,'full').*wwm2)>=1,c3,'same')<1)); % particles mask
%       pmask=wp((c2x-1)/2+1:end-(c2x-1)/2,(c2y-1)/2+1:end-(c2y-1)/2).*(1-ws);
% 
%       img=pmask.*img;
% 
%       chiimg=1+(-2*conv2l(img,smip.^2,'full')+conv2l(img.^2,smip,'full'))./ww3;
%       [Np,spx,spy]=findpeakstr(chiimg,kern,mk,GoodBallsCutOff,MinSepFrac*D);
%       ri=zerofill(img,NNx-Nx,NNy-Ny);
%       [spx,spy]=refinepos(ri,spx,spy,rect,D,w,dD,edgmask);
%       Dw=fminsearch('chiDw',[D w],optimset('TolX',1e-2,'TolFun',10),ri,spx,spy,rect,dD,Np,edgmask,hs,0);
%       D=Dw(1)
%       w=Dw(2)
% 
%       [dD,xxs,yys,smip,ww3,mk,NNx,NNy,sx,sy,rect,hs,hp,c1,wwm1,c2,wwm2,c3,edgmask]=DwdependentVar(D,w,Nx,Ny,static,refine);
% 
%       chiimg=1+(-2*conv2l(img,smip.^2,'full')+conv2l(img.^2,smip,'full'))./ww3;
%       [Np,spx,spy]=findpeakstr(chiimg,kern,mk,GoodBallsCutOff,MinSepFrac*D);
%       ri=zerofill(img,NNx-Nx,NNy-Ny);
%       [spx,spy]=refinepos(ri,spx,spy,rect,D,w,dD,edgmask);
%       Dw=fminsearch('chiDw',[D w],optimset('TolX',1e-4,'TolFun',10),ri,spx,spy,rect,dD,Np,edgmask,hs,0);
%       D=Dw(1)
%       w=Dw(2)
% 
%       fprintf(1,'Saving diameter and width...\n');
%       save(Dwfn,'D','w');
% 
%       [dD,xxs,yys,smip,ww3,mk,NNx,NNy,sx,sy,rect,hs,hp,c1,wwm1,c2,wwm2,c3,edgmask]=DwdependentVar(D,w,Nx,Ny,static,refine);
%     end
%   end
%   fprintf(1,'Done.\n');
% 
%   [c1x c1y]=size(c1);
%   [c2x c2y]=size(c2);
% 
%   t0=clock;
% 
%   sdkmk=0;
%   sbtmk=0;
%   dkbg=0;
%   btbg=0;
%   for nf=[resume:resume_end];
%     raw=readimgarr(fnimg,Nx+clipx,Ny+clipy,frames(nf));
%     img=clip((bg1-raw)./bg2.*(1-static),mncout,mxcout);
%     img=img(clipx+1:end,clipy+1:end);% clip crud at side
%     if(min(img(:))>.2)
%       mni=min(img(:));
%       img=(img-mni)/(1-mni);
%     end
%     ws=(1-(conv2l((conv2l(img>ithresh,c1,'full').*wwm1)>=1,c1,'same')<1)); % wall mask
%     ws=ws((c1x-1)/2+1:end-(c1x-1)/2,(c1y-1)/2+1:end-(c1y-1)/2);
%     wp=(1-(conv2l((conv2l((img>ithresh)-ws,c2,'full').*wwm2)>=1,c3,'same')<1)); % particles mask
%     pmask=wp((c2x-1)/2+1:end-(c2x-1)/2,(c2y-1)/2+1:end-(c2y-1)/2).*(1-ws);
%     img=pmask.*img;
%     chiimg=1+(-2*conv2l(img,smip.^2,'full')+conv2l(img.^2,smip,'full'))./ww3;
%     [Npf(nf),spx,spy]=findpeakstr(chiimg,kern,mk,GoodBallsCutOff,MinSepFrac*D);
%     if(length(spx)>0)
%       if (refine==1)
%         ri=zerofill(img,NNx-Nx,NNy-Ny);
%         ntimes=0;
%         [spx,spy]=refinepos(ri,spx,spy,rect,D,w,dD,edgmask);
%         ci=calcimg(spx,spy,NNx,NNy,rect,Npf(nf),D,w,dD,1);
%         sc=sum(ci(:).^2)./sum(ci(:).*ri(:));
%         box=(sc*ri-ci);
%         di3=-box.*edgmask;
%         chi(nf)=sum(di3(:).^2);
%         while(((chi(nf)/Npf(nf)>.6) && (ntimes<3)) )
%           ntimes=ntimes+1;
%           [spx,spy]=refinepos(ri,spx,spy,rect,D,w,dD,edgmask);
%           ci=calcimg(spx,spy,NNx,NNy,rect,Npf(nf),D,w,dD,1);
%           sc=sum(ci(:).^2)./sum(ci(:).*ri(:));
%           box=(sc*ri-ci);
%           di3=-box.*edgmask;
%           chi(nf)=sum(di3(:).^2);
%         end
%       end
% 
%       px(1:Npf(nf),nf)=spx;
%       py(1:Npf(nf),nf)=spy;
% 
%       if (plotit==1 & rem(nf-1,nplot+1)==0)
%         %set(hs,'Cdata',zerofill(img,NNx-Nx,NNy-Ny));
%         %    set(hs,'Cdata',zerofill(img,NNx-Nx,NNy-Ny));
%         if (refine==1)
%           set(hs,'Cdata',[ri di3]);
%           caxis([-.25 .25]);
%         else
%           set(hs,'Cdata',zerofill(img,NNx-Nx,NNy-Ny));
%           caxis([-inf inf]);
%         end
%         set(hp,'xdata',py(1:Npf(nf),nf),'Ydata',px(1:Npf(nf),nf));
%         drawnow;
%         if(pauseit==1)
%           pause;
%         end
%       end
%       dkmk=calcimg(spx-(NNx-Nx)/2+1,spy-(NNy-Ny)/2+1,Nx+1,Ny+1,[1 Nx+1 1 Ny+1],Npf(nf),D*(1-.08),10,dD,1);
%       btmk=1-calcimg(spx-(NNx-Nx)/2+1,spy-(NNy-Ny)/2+1,Nx+1,Ny+1,[1 Nx+1 1 Ny+1],Npf(nf),D*(1+.08),10,dD,1);
%       sdkmk=sdkmk+dkmk;
%       sbtmk=sbtmk+btmk;
%       dkbg=dkbg+raw.*dkmk;
%       btbg=btbg+raw.*btmk;
%       figure(2);
%       simage(btbg./(sbtmk+eps));
%       figure(3);
%       simage(dkbg./(sdkmk+eps));
%       drawnow;
%       figure(1);
%     end
%     done(nf)=1;
%     if(rem(nf-1,100)==0)
%       fprintf('saving progress...');
%       tpx=px;
%       tpy=py;
%       NBm=max(Npf(:));
%       px(NBm+1:end,:)=0;
%       py(NBm+1:end,:)=0;
%       save(trfn,'px','py','Npf','pbox','chi','done');
%       px=tpx;
%       py=tpy;
%       fprintf('Done\n');
%     end
% 
%     disp([nf/10 etime(clock,t0)/60 etime(clock,t0)*((Nf-resume+1)/((nf-resume+1))-1)/60 Npf(nf) str2num(id)]);
%     diary off;
%     diary([fn '.diary']);
%   end
% 
%   if (saveit==1)
%     if(exist(trfn))
%       delete(trfn);
%     end
%     NBm=max(Npf(:));
%     px(NBm+1:end,:)=0;
%     py(NBm+1:end,:)=0;
%     sfn=sprintf('%s%s_%s.mat','track',rev,fn);
%     save(sfn,'px','py','Npf','pbox','chi','done','info');
%   end
% 
% catch
%   disp('error capture:');
%   rethrow(lasterror);
%   diary off;
%   diary([fn '.diary']);
% end
% 
% 
% 
% 
% 
