%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';
  fn='Vert2D_Run54_070717_192158_0080';
 % fn='Vert2D_Run54_070717_191150_0050';
 % fn='Vert2D_Run54_070717_185820_0010'
  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','static','Dw','D','DD','w','ww','px','py','Npf','done','chi','pbox','tbox'};  % 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=7.6;                    % first guess for D
  w=1;                    % first guess for w
  MaxBalls=140;            % Maximum number of balls to find
  NB=136;
  GoodBallsCutOff=5;       % Cutoff peak height for ball (not D independent)
  FracPart=.1;            % 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=5;       % width of mask to lookfor second peak
  mnmxfrac=.01;     % Nbg*mnmxfrac lowest and Nbg*(1-mnmxfrac) highest used to calculate mn and mx
  staticfrac=.01;   % 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
  pbox=zeros(6,Nf);        % particle bounding box
  tbox=zeros(6,Nf);        % temp particle bounding box
  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)+1;                                               % position offset relative to raw image
      oy=fix((dD-clipy)/2)+1;
      [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];
     
      x1=max([1 find(diff(mean(static,1)<.995)==1,1)]);
      x2=find(diff(mean(static,1)<.995)==-1,1,'last');
      yy1=find(diff(mean(static,2)<.9)==1,1);
      yy2=find(diff(mean(static,2)<.9)==-1,1,'last');
      pbox(5:6,:)=repmat([x1;x2],1,Nf);
      sNy=yy2-yy1+1;

      if(plotit==1)  % plot progress
        hs=simage(zeros(3*sNy,Nx));                
        hold on
        hp=plot(1:2,1:2,'k.');
        hold off;
      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;
        ri=img;
        y1=find(diff(mean(ri,2)<.95)==1,1);
        y2=find(diff(mean(ri,2)<.95)==-1,1,'last');
        pbox(1:2,nf)=y1;
        pbox(3:4,nf)=y2;
        bmask=zeros(size(img));
        tt=img(y1:y2,x1:x2);
        img=bmask;
        img(y1:y2,x1:x2)=tt;
        bmask(y1+4:y2-3,x1+2:x2-2)=1;
        
        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  
          %pause;
          [Npf(nf),spx,spy]=findpeaks(1./chiimg,mk,GoodBallsCutOff,MinSepFrac*D);      % find all maxima > CutOff with separation > Frac*D
          if(Npf(nf)>NB)
            Npf(nf)=NB;
            spx=spx(1:NB);
            spy=spy(1:NB);
          end
          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.
        di=bmask.*(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=bmask.*(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=bmask.*(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=bmask.*(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=bmask.*(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=bmask.*(ci-ri);

        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);

        box=ri-ci;
        D2=fix(D/2);
        ix=pbox(5,nf)+D2:pbox(6,nf)-D2;
        iytop=fix(mean(pbox(1:2,nf)))+(-D2:D2);
        [mm btop]=min(diff(box(iytop,ix)));
        iybot=fix(mean(pbox(3:4,nf)))+(-D2:D2);
        [mm bbot]=max(diff(box(iybot,ix)));
        hold off;
        tbox(:,nf)=pbox(:,nf);
        pbox(1:2,nf)=polyval(polyfit(ix,iytop(btop),1),[1 Nx])';
        pbox(3:4,nf)=polyval(polyfit(ix,iybot(bbot),1),[1 Nx])';
        

        if (plotit==1 & rem(nf-1,nplot+1)==0)
%          set(hs,'Cdata',[bgmx./(nmx+eps) raw bgmn./(nmn+eps)]);
          set(hs,'Cdata',[di(yy1:yy2,:); ri(yy1:yy2,:); di2(yy1:yy2,:)]);
          caxis([-.4 .4]);
          set(hp,'xdata',spy-oy,'Ydata',spx-ox+sNy-yy1+1);
          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
  
