%% Track

ver=ver();  % set version

%% User Input
idir='ImageRepository';
fns='*.JPG'; % filenames to process

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

D=13;       % Particle diameter
w=1;       % Particle edge width
CutOff=300;  % Particle detection cutoff 
MT=260;     % Max tetramer
MPT=4;      % Max particle per tetramer

refine=false;

recalcbg=false;  % Force background recalc
cut=200;        % Background cutoff

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

%% Calculated Parameters

% File information
dr=dir([idir '/' fns]);
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'];

% Background and time and box
prfn=[base '_pre.mat'];
if(~exist(prfn,'file') || recalcbg)
	trackprep;
else
	load(prfn);
end
	colormap(jet(256));

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

% background
mx=mx(box(1):box(2),box(3):box(4),1)/3;
mn=mn(box(1):box(2),box(3):box(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;
frames=sf:ef;
px=zeros(MT,MPT,Nf);
py=zeros(MT,MPT,Nf);
Npf=zeros(1,Nf);
sNpf=zeros(MT,Nf);
[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;
ppx=zeros(MT,MPT,Nf);
ppy=zeros(MT,MPT,Nf);

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

%% Main Loop
Ds=D;
ws=w;
tic;
for nf=1:Nf
	% load frame
	nff=frames(nf);
	im=double(imread([idir '/' dr(nff).name]));
  
	% roi
	raw=im(box(1):box(2),box(3):box(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,1,CutOff,30/4);
	Npf(nf)=min(Npf(nf),MT);
	spx=4*spx(1:Npf(nf))-2;
	spy=4*spy(1:Npf(nf))-2;
	px(1:Npf(nf),1,nf)=spx;
	py(1:Npf(nf),1,nf)=spy;
	
	for np=1:Npf(nf);
		D=Ds;
		w=ws;
		ims=img(px(np,1,nf)+(-22:22),py(np,1,nf)+(-22:22));
		[sNp spx spy]=findpeaks(1./chiimg(ims-sci,ip,[],sw3),mks,4,D);
		if(sNp>3)
			sNp=3;
			spx=spx(1:3);
			spy=spy(1:3);
		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); colorbar;
	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],varsavlist{:});
		fprintf('Done.\n');
	end
end
toc;

fprintf('Saving...');
save(trfn,varsavlist{:});
fprintf('Done.\n');

