%% Track

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

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

D=11;       % Particle diameter
w=1;       % Particle edge width
CutOff=10;  % Particle detection cutoff 
MB=1100;    % Max particles

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

varsavlist={'px','py','Npf','ppx','ppy',...
	          'nns','sclist','mbr','D','w',...
						'imsc','frames','base'};

%% Calculated Parameters

% File information
dr=dir(fns);
info=imfinfo(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)
	trackprep000;
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>
[bNx bNy]=size(mx);

% 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

%% Define

Nf=ef-sf+1;
frames=sf:ef;
px=zeros(MB,Nf);
py=zeros(MB,Nf);
Npf=zeros(1,Nf);
imsc=zeros(1,Nf);
mbr=zeros(1,Nf);
[junk w3]=chiimg(mx,ip);
mk=conv2(0*mx+1,ip,'full')./sum(ip(:))>.5;
nns=zeros(Nf,13);
sclist=zeros(fix(MB/4),Nf);
ppx=zeros(fix(MB/4),Nf);
ppy=zeros(fix(MB/4),Nf);

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

%% Main Loop

for nf=1:Nf
	nff=frames(nf);
	im=double(imread(dr(nff).name));
	raw=im(box(1):box(2),box(3):box(4),1);
	br=(mx-raw)<60;
	mbr(nf)=mean(br(:));
	imsc(nf)=mean(raw(br))/mean(mx(br));
	if(imsc(nf)>.85 && mbr(nf)>.85)
		img=1-raw./mx/imsc(nf);
		[nn bb]=hist(img(:),bb);
		[junk2 p1]=max(nn(1:bbm));
		[junk1 p2]=max(nn(bbm+1:end));
		img=clip((img-bb(p1))/(bb(p2+bbm-1)-bb(p1)),-.1,1.1);
		[Npf(nf) spx spy]=findpeaks(1./chiimg(img,ip,[],w3),mk,CutOff,.75*D);
		Npf(nf)=min(Npf(nf),MB);
		spx=spx(1:Npf(nf));
		spy=spy(1:Npf(nf));
		px(1:Npf(nf),nf)=spx;
		py(1:Npf(nf),nf)=spy;
		Xn=repmat(spx+i*spy,1,Npf(nf)); %#ok<IJCL>
		dd=abs(Xn-Xn.');
		clist=zeros(1,Npf(nf));
		for np=1:Npf(nf)
			ii=find(dd(np,:)<50);
			[junk pp]=min(abs(Xn(ii)-mean(Xn(ii))));
			clist(np)=ii(pp);
		end
		clist=unique(clist);
		Ntet=min(length(clist),256);
		ppx(1:Ntet,nf)=spx(clist(1:Ntet));
		ppy(1:Ntet,nf)=spy(clist(1:Ntet));
		simage(img);
		hold on; 
		plot(spy-ss,spx-ss,'r.');
		plot(ppy(1:Ntet,nf)-ss,ppx(1:Ntet,nf)-ss,'w.');
		axis('image');
		hold off;
		drawnow;
		fprintf('%d %d\n',nf,Npf(nf));
	else
		Npf(nf)=-99;
	end
	if(mod(nf,100)==0)
		fprintf('Saving...');
		save(['tmp_' trfn],varsavlist{:});
		fprintf('Done.\n');
	end
end

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