%% Track

%% User Input
idir='ImageRepository';  % relative directory for images
fns='*.JPG'; % filenames to process

pt=5;    % time to pause to check for correct preprocessing

bd=50;      % boundary of cell to exclude
pac=8;      % paraxial correction to bd
            % (Nx-2bd) must be divisible by 4
            % (Ny-2bd+pac) must be divisible by 4
Nb=200;     % number of randomly choosen frames to get background
froff=100;  % frame offset to get background (nb: not longer used)

sf=1;        % start frame
ef=[];       % end frame []==all

D=13;       % Particle diameter
w=1;        % Particle edge width

CutOffTet=300;  % Whole tetramer detection cutoff
CutOffPart=4;   % Individual particle detection cutoff
MT=256;      % Max # of tetramer
MPT=4;       % Max # particle per tetramer

refine=false;  % does not work

restart=true;  % true trys to restart from tmp_track file

recalcbg=false;  % Force background recalc

% List of variables to save

tracksavlist={'px','py','Npf','ppx','ppy','sNpf',...
    'D','w',...
    'frames','base','ver','info'};

timesavlist={'time','segt','seg','ledt','led','segraw','ledraw','base'};
clksavlist={'segii','ledii','seglut','segx','segdx','segy','segdy','rx','ry'};

ver=ver();  % set version

overide('track',[],'_user'); % override defaults
%% Calculated Parameters

% File information
fprintf('Processing Image Filenames...');
dr=dir([idir '/' fns]);
fprintf('Done.\n');
Nff=length(dr);
info=imfinfo([idir '/' dr(1).name]);
base=info.DateTime([3 4 6 7 9 10 12 13 15 16 18 19]);
trfn=['track_' base '.mat'];
bgfn=['bg_' base '.mat'];
clkfn=['clk_' base '.mat'];
timefn=['time_' base '.mat'];

% Background and box
if(~exist(bgfn,'file') || recalcbg)
		fprintf('Calculating Background...');
		rfr=fix(Nff*rand(1,2*Nb))+1;
		rfr=unique(rfr);
		[junk ii]=sort(rand(1,length(rfr)));
		rfr=sort(rfr(ii(1:Nb)));
		aimg=0;
		mx=0;
		mn=9e99;
		Nc=length(sprintf('%d',Nb));
		fmt=sprintf('[%%0%dd/%%0%dd]',Nc,Nc);
		fprintf(1,fmt,0,Nb);
		for nb=1:Nb
			nff=rfr(nb);
			im=double(imread([idir '/' dr(nff).name]));
      aimg=aimg+im(:,:,1);    
      mx=max(mx,im(:,:,1));    
      mn=min(mn,im(:,:,1));
			h=image(uint8(im));
			set(h,'cdatamapping','direct');
			axis('image');
			drawnow;
			backspace(1,fmt,nb,Nb);
			fprintf(1,fmt,nb,Nb);
		end
		aimg=3*aimg/Nb;  % 3 is for compatibility
		mx=mx*3;
		mn=mn*3;
		colormap(jet(256));
		fprintf('...Done.\n');
		fprintf('Saving Background...');
    save(bgfn,'aimg','mx','mn','Nff','info','base');
		fprintf('Done.\n');
else
		fprintf('Loading Background...');
    load(bgfn);
		fprintf('Done.\n');
end

info.xls=parsexls();  % parse xls file

if(~exist('bbox','var'))
    bbox=findbox(aimg);
    bbox=bbox+[bd -bd +bd -bd+pac];  % exclude boundary
    
    simage([mx(bbox(1):bbox(2),bbox(3):bbox(4)) mn(bbox(1):bbox(2),bbox(3):bbox(4))]);
    drawnow;
    
    fprintf('Use ctrl-c to stop, if background or box are not correct.\n');
    countdown(pt);
    save(bgfn,'bbox','aimg','mx','mn','Nff','info','base','bd','pac');
end

if(~exist(clkfn,'file'))
	% LED 7-SEG clock
	[segii ledii seglut segx segdx segy segdy rx ry Nl]=...
		getledclkparm(aimg,bbox,bd,[idir '/' dr(froff).name]);
  if(Nl~=8)
		[segii ledii seglut segx segdx segy segdy rx ry Nl]=...
			getledclkparm(mx,bbox,bd,[idir '/' dr(froff).name],[],700);
	end
	% check LED 7-seg data
	simage(aimg(reshape(ledii,11,11*8)));
	drawnow;
	fprintf('Use ctrl-c to stop, if LEDs are not correct.\n');
	countdown(pt);
	
	simage(aimg(:,:,1));
	axis([segy+1.2*segdy/2*[-1 1],segx+segdx/2*[-1 1]]);
	hold on;
	plot(segy+ry,segx+rx,'*');
	hold off;
	drawnow;
	fprintf('Use ctrl-c to stop, if 7-seg are not correct.\n');
	countdown(pt);
  save(clkfn,clksavlist{:});
else
	load(clkfn);
end

% defaults
if(isempty(ef))
    ef=Nff;
end

% background
mx=mx(bbox(1):bbox(2),bbox(3):bbox(4),1)/3;
mn=mn(bbox(1):bbox(2),bbox(3):bbox(4),1)/3; %#ok<NASGU>
[Nx Ny]=size(mx);

% smoth background
mx1=mx;
mx=conv2(mx,ones(5,5),'same')./conv2(ones(size(mx)),ones(5,5),'same');

% Ideal Particle
ss=fix(D/2)+1;  % Ideal particle box size

[xx yy]=ndgrid(-ss:ss,-ss:ss);
rr=xx+i*yy; %#ok<IJCL>
ip=ipf(abs(rr),D,w);  % Ideal particle

% Tetramer mask
sst=fix(D/4)+1;  %
[xx yy]=ndgrid(-sst:sst,-sst:sst);
rr=xx+i*yy; %#ok<IJCL>
tmk=double(abs(rr)<sst+1);  % Tetramer mask

%% create local grid;
[cxy over]=pgrid(23,23,45,45,[1 45 1 45],1,3*D,0);
sci=ipf(cxy,D,w);


%% Define

Nf=ef-sf+1;        % # of files
frames=sf:ef;      % list of frames

px=zeros(MT,MPT,Nf);  % crude particle positions
py=zeros(MT,MPT,Nf);
ppx=zeros(MT,MPT,Nf); % refined particle positions
ppy=zeros(MT,MPT,Nf);

Npf=zeros(1,Nf);      % Number of tetramers
sNpf=zeros(MT,Nf);    % Number of particles/tetramer

% fitting defines
[junk w3]=chiimg(mx,ip); %#ok<*ASGLU>
[junk sw3]=chiimg(sci,ip);
mk=conv2(0*mx+1,ip,'full')./sum(ip(:))>.5;
mks=conv2(0*sci+1,ip,'full')./sum(ip(:))>.5;

% LED clock defines
segraw=zeros(7,4,Nf);
ledraw=zeros(8,Nf);

% misc defines
bb=-.5:.01:1.5;
[mn bbm]=min(abs(bb-.5));

rnd=1e-8*randn(Nx/4,Ny/4); % add to chi so no pix are equal

%% Main Loop
Ds=D;
ws=w;
tic;
if (restart && exist(['tmp_' trfn],'file'))
	fprintf('Loading restart info...');
	load(['tmp_' trfn]);
	snf=find(Npf==0,1);
	fprintf('Done.\n');
else
	snf=sf;
end
for nf=snf:Nf
    % load frame
    nff=frames(nf);
    im=double(imread([idir '/' dr(nff).name]));
    
    % raw time calc
    segraw(:,:,nf)=squeeze(sum(im(segii)));
    ledraw(:,nf)=sum(im(ledii));
    
    % roi
    raw=im(bbox(1):bbox(2),bbox(3):bbox(4),1);
    
    % normalize frame
    img=raw./mx;
    img=histnorm(img,bb,bbm);
    img=clip(img,-.15,1.15,0,1);
    
    % peak at tetramer position
    chi=permute(reshape(img.*(img>.2),4,Nx/4,4,Ny/4),[2 4 1 3]);
    chi=conv2(sum(sum(chi,4),3),tmk,'same');
    
    % find tetramers
    [Npf(nf) spx spy]=findpeaks(chi+rnd,1,CutOffTet,30/4);
    Npf(nf)=min(Npf(nf),MT);
    spx=4*spx(1:Npf(nf))-1.5;
    spy=4*spy(1:Npf(nf))-1.5;
    px(1:Npf(nf),1,nf)=spx;
    py(1:Npf(nf),1,nf)=spy;
    
    for np=1:Npf(nf);
        D=Ds;
        w=ws;
        lx=clip(ceil(px(np,1,nf))+(-22:22),1,Nx);
        ly=clip(ceil(py(np,1,nf))+(-22:22),1,Ny);
        ims=img(lx,ly);
        sNp=0;
        fc=1;
        while(sNp<3 && fc>.6);
            [sNp spx spy]=findpeaks(1./chiimg(ims-sci,ip,[],sw3),mks,CutOffPart*fc,D);
            if(sNp>3)
                sNp=3;
                spx=spx(1:3);
                spy=spy(1:3);
            end
            fc=fc*.95;
            if(sNp<3)
                fprintf('.');
            end
        end
        spx=[23;spx-ss];
        spy=[23;spy-ss];
        sNp=sNp+1;
        sNpf(np,nf)=sNp;
        px(np,1:sNp,nf)=spx+px(np,1,nf)-23;
        py(np,1:sNp,nf)=spy+py(np,1,nf)-23;
        [cxy over]=pgrid(spx,spy,45,45,[1 45 1 45],sNp,3*D,0);
        ci=ipf(cxy,D,w);
        ims=ims.*ipf(cxy,D*1.2,w);
        di=ci-ims;
        [dpx dpy]=cidp2(cxy,over,di,sNp,D,w);
        spx=spx+dpx;
        spy=spy+dpy;
        if(refine)
            [cxy over]=pgrid(spx,spy,45,45,[1 45 1 45],sNp,3*D,0);
            ci=ipf(cxy,D,w);
            di=ci-ims;
            [delD delw]=cidDw(abs(cxy),di,D,w);
            D=D+delD;
            w=w+delw;
            ci=ipf(cxy,D,w);
            di=ci-ims;
            [dpx dpy]=cidp2(cxy,over,di,sNp,D,w);
            spx=spx+dpx;
            spy=spy+dpy;
        end
        ppx(np,1:sNp,nf)=spx+px(np,1,nf)-23;
        ppy(np,1:sNp,nf)=spy+py(np,1,nf)-23;
    end
    
    % user feedback
    simage(img);
    hold on;
    plot(flat(ppy(1:Npf(nf),:,nf)),flat(ppx(1:Npf(nf),:,nf)),'w.');
    axis('image');
    hold off;
    drawnow;
    fprintf('%d %d %d\n',nf,Npf(nf),sum(sNpf(:,nf)));
    if(mod(nf,100)==0)
        hold off
        plot(flat(ppx(1:256,:,nf-99:nf)-ppx(1:256,ones(1,4),nf-99:nf)),flat(ppy(1:256,:,nf-99:nf)-ppy(1:256,ones(1,4),nf-99:nf)),'.')
        hold on
        plot(flat(px(1:256,:,nf-99:nf)-px(1:256,ones(1,4),nf-99:nf)),flat(py(1:256,:,nf-99:nf)-py(1:256,ones(1,4),nf-99:nf)),'r.')
        axis('equal');
        hold off;
        fprintf('Saving...');
        save(['tmp_' trfn],tracksavlist{:});
        fprintf('Done.\n');
    end
end
toc;

%% Save results

fprintf('Saving final track infomation...');
save(trfn,tracksavlist{:});
fprintf('Done.\n');

fprintf('Processing raw time information...');
[time segt seg ledt led]=procrawclk(segraw,ledraw,seglut);
fprintf('Done.\n');

fprintf('Saving time infomation...');
save(timefn,timesavlist{:});
fprintf('Done.\n');

