function [GRP,PARMS]=gram_detect(CFAR,jdim_proc,cfar_thresh)
%
% Call [GRP,PARMS]=gram_detect(CFAR,jdim_proc,cfar_thresh);
%
%
% COPYRIGHT (c) 2008
% 
% Applied Research Laboratories, The University of Texas at Austin
% (ARL:UT)
% All rights reserved.
% 
% This material may be reproduced by or for the U.S. Government
% only. This notice must appear in all copies of this file and its
% derivatives. Further distribution outside the U.S. Government is
% prohibited without the explicit written authorization from the
% Director of the Environmental Sciences Laboratory, ARL:UT.
% 
% THIS SOFTWARE IS PROVIDED "AS IS" WITHOUT WARRANTY OF ANY KIND,
% EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO THE
% IMPLIED WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
% PURPOSE, OR NONINFRINGEMENT.
% 
% IN NO EVENT WILL ARL:UT BE LIABLE FOR ANY DAMAGES WHATSOEVER ARISING
% OUT OF THE USE OR INABILITY TO USE THIS SOFTWARE, EVEN IF ARL:UT HAS
% BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGES.
if nargin<3, cfar_thresh=5, end
%
view_global;
%
% Additional parameters:
npt_min=1; cfar_delta=2;
%
jdim1=jdim_proc(1); jdim2=jdim_proc(2);
%
ndim_loop=0; nloop=1;
for jdim=1:CFAR.ndim,
    if jdim~=jdim1 & jdim~=jdim2, %% & CFAR.dims(jdim)>1,
        ndim_loop=ndim_loop + 1;
        jdim_loop(ndim_loop)=jdim;
        dims_loop(ndim_loop)=CFAR.dims(jdim);
        nloop=nloop * CFAR.dims(jdim);
        indx_dims_full{ndim_loop}=[1:CFAR.dims(jdim)];
    end
end
%
if ndim_loop==0,
    disp('CFAR has only the two dimensions');
    ndim_loop=1; jdim_loop=[]; dims_loop=[]; nloop=1; % Case where a single image exists;
end;     
%
ntime=CFAR.dims(jdim1); nbrg=CFAR.dims(jdim2);
t_ax=CFAR.scales{jdim1}; brg_ax=CFAR.scales{jdim2};
iimid=[1:nbrg]; iilo=[1,1:nbrg-1]; iihi=[2:nbrg,nbrg];
%
dims_loop_over=CFAR.dims(jdim_loop); 
idx=cell(2+ndim_loop,1);      % idx is cell array of indices into CFAR.data;
idx{jdim1}=[1:CFAR.dims(jdim1)]; idx{jdim2}=[1:CFAR.dims(jdim2)];   % Fill in the two image dims;
%
npt_max=10*ntime; npt_max0=0;
ng_max=100; ng_max0=0; ng_unique0=0; blob_parms=zeros(ng_max,6);
GRP.ndim=2 + ndim_loop;
GRP.dims=[ntime,ng_max,dims_loop];
GRP.labels={CFAR.labels{jdim1}, 'Group #', CFAR.labels{jdim_loop}};
GRP.scales={CFAR.scales{jdim1}, [1:ng_max], CFAR.scales{jdim_loop}};
GRP.data=NaN*ones(GRP.dims);
GRP.dlab=CFAR.labels{jdim2};
if ndim_loop>1, GRP.n_group=zeros(dims_loop); 
elseif ndim_loop==1, GRP.n_group=zeros(ndim_loop,1);
elseif ndim_loop==0, GRP.n_group=0;
end
%
PARMS.ndim=2 + ndim_loop;
PARMS.dims=[ng_max,6,dims_loop];
PARMS.labels={'Group No', 'Parm #', CFAR.labels{jdim_loop}};
PARMS.scales={[1:ng_max] [1:6], CFAR.scales{jdim_loop}};
PARMS.dlab='Slope/t0/CFmax/Delta';
%
slope_fac=1000; t_fac=60;      % Convert slope, time to units of choice;
fit_fac=3; 
for jloop=1:nloop,
    if ndim_loop>0, indx_dims=indx_full_calc(dims_loop_over,jloop); end
    for j=1:ndim_loop, 
        idx{jdim_loop(j)}=indx_dims(j); 
        idx_pad{j}=indx_dims(j);
    end
%    disp(idx)
    data=squeeze(CFAR.data(idx{:}));
    if jdim1<jdim2, data=data.'; end
%    
% Now data is nbrg by ntime:
	iimax1=find(data(iimid,:)>cfar_thresh & data(iimid,:)>data(iilo,:) & data(iimid,:)>data(iihi,:));
	indx_full=indx_full_calc([nbrg,ntime],iimax1);
	nmax=length(iimax1); brg_max=zeros(nmax,1); t_max=brg_max; jt_max=brg_max;
	for j=1:nmax,
		iib=[indx_full(j,1)-1:indx_full(j,1)+1];
        data_vals=data(iib,indx_full(j,2));
		abc=polyfit(brg_ax(iib)',data_vals,2);
		brg_max(j)=-abc(2)/(2*abc(1));
        cfar_max(j)=abc(1)*brg_max(j)^2 + abc(2)*brg_max(j) + abc(3);
        t_max(j)=t_ax(indx_full(j,2));
        jt_max(j)=indx_full(j,2);
	end
%
	dd_thresh0=16.5;
	ngroup_max=100; nmember_max=100;
	nmember=zeros(1,ngroup_max); kmember=zeros(nmember_max,ngroup_max);
	ngroup=1; nmember(ngroup)=1; kmember(1,ngroup)=1;
	iidiag=0;
	for j=2:nmax
        jbrg=indx_full(j,1); jtime=indx_full(j,2);
        if iidiag, disp(['Testing ',num2str([jbrg,jtime])]); end
        dd_min=1e20;
        for jgroup=1:ngroup
% Test last member added to each group:
            jg=kmember(nmember(jgroup),jgroup); 
            jbrg_test=indx_full(jg,1); jtime_test=indx_full(jg,2);
            dd=(jbrg-jbrg_test).^2 + (jtime - jtime_test).^2;
            if iidiag, disp([jbrg,jbrg_test,jtime,jtime_test,dd,dd_min]); end
            if dd<dd_min, dd_min=dd; jg_min=jgroup; end
        end
        if dd_min>dd_thresh0,
            if iidiag, disp(['Group added: dd_min = ',num2str(dd_min)]); end
            ngroup=ngroup + 1;
            nmember(ngroup)=1;
            kmember(nmember(ngroup),ngroup)=j;
        else
            if iidiag, disp(['Member added to group ',int2str(jg_min)]); end
            nmember(jg_min)=nmember(jg_min) + 1;
            kmember(nmember(jg_min),jg_min)=j;
        end
	end
%
	if ngroup==0, disp(['Number of groups found to be zero']); end
%
	if iidiag,
		for jg=1:ngroup,
            disp(['Group ',int2str(jg),': Number of members = ',int2str(nmember(jg))]);
		end
	end
%	figure(101); clf; 
cols='rbgmck'; ncol=length(cols); zoom on;
%
% Edit out groups consisting of too few members or with too steep or shallow slopes:
	ng_ok=0; ng_unique=0;
	slp_min=1; slp_max=20;
	kmemb_ok=zeros(1,ngroup_max); kg_unique=zeros(1,ngroup_max);
    Pvals=zeros(ngroup_max,2);
    for jg=1:ngroup,
        iiok=1;
        if nmember(jg)<=npt_min, iiok=0; end
        if iiok,
            km=kmember(1:nmember(jg),jg);
            tvals=t_fac*t_max(km); brgvals=brg_max(km);
            [cf_max,ii_max]=max(cfar_max(km)); time_max=tvals(ii_max); bearing_max=brgvals(ii_max);
            [P,S]=polyfit(tvals,brgvals,1);
            t0=-P(2)/P(1);
            slp=P(1)*slope_fac;      % convert to units of choice;
            if abs(slp)<slp_min | abs(slp)>slp_max, iiok=0; end
% Require one point to exceed cfar_thresh+cfar_delta;
%            if cf_max<cfar_thresh+cfar_delta, iiok=0; end       
%            [y0,delta0]=polyval(P,t0,S); delta_t0=3*abs(delta0/P(1));
            disp(['New group: ',num2str([t0,slp])]);
        end
        if iiok
            if iidiag, disp(['Group added: dd_min = ',num2str(dd_min)]); end
            ng_ok=ng_ok + 1; 
            if ng_ok>ng_max, disp('Increas ng_max in gram_detect!!'); error('Exiting'); end
            kmemb_ok(ng_ok)=ng_ok;
%            [y,delta]=polyval(P,tvals,S);
            Pvals(ng_ok,:)=P; Svals(ng_ok)=S;
            blob_parms(ng_ok,1)=slp;
            blob_parms(ng_ok,2)=t0;
            blob_parms(ng_ok,3)=cf_max;
            blob_parms(ng_ok,4)=time_max;
            blob_parms(ng_ok,5)=bearing_max;
            blob_parms(ng_ok,6)=1;
%           
            t0_thresh=5; slp_thresh=3;
            indx_grp={[], [], idx_pad{:}};
            for jgroup=ng_ok-1:-1:1,
% Test this group with all other groups to see if slope and intercept change too much:
                jg2=kmemb_ok(jgroup);
% Extrapolate forward from current group to max of group jg2:
                [brg_extrap1,delta1]=polyval(P,blob_parms(jg2,4),S);
% Extrapolate backward from group jg2 to max of current group:
                [brg_extrap2,delta2]=polyval(Pvals(jg2,:),time_max,Svals(jg2));
                fit1=abs(brg_extrap1 - blob_parms(jg2,5)) / (fit_fac*delta1);
                fit2=abs(brg_extrap2 - bearing_max) / (fit_fac*delta2);
                disp(['Extrap1: ',num2str([blob_parms(jg2,4),brg_extrap1,blob_parms(jg2,5),delta1,fit1])]);
                disp(['Extrap2: ',num2str([time_max,brg_extrap2,bearing_max,delta2,fit2])]);
                if fit1 + fit2 < 2.5 & slp*blob_parms(jg2,1)>0,
                    disp(['Association between ',int2str([jg,jg2])]);
                    iiok=0;
                    break;
                else
                    disp(['Association rejected.']);
				end
            end
            if iiok==0,
                jg_unique=kg_unique(jg2);
                kg_unique(ng_ok)=jg_unique;
				indx_grp={jt_max(km), jg_unique, idx_pad{:}};
				GRP.data(indx_grp{:})=brgvals;
                disp(['Joining group with unique group.',int2str(jg_unique)]);
            else
                disp(['Entering unique group.']);
                ng_unique=ng_unique + 1;
                kg_unique(ng_ok)=ng_unique;
				indx_grp={jt_max(km), ng_unique, idx_pad{:}};
				GRP.data(indx_grp{:})=brgvals;
             end
        end
	end
%
    ng_max0=max(ng_ok,ng_max0),
    ng_unique0=max(ng_unique,ng_unique0);
%
    indx_parms={[1:ng_ok], [1:6], idx_pad{:}};
	PARMS.data(indx_parms{:})=blob_parms(1:ng_ok,:);
	%
	figure(102); clf;
	V.ii_img=1; vdims=[jdim1,jdim2,jdim_loop];
	view5(CFAR,vdims,0,0); hold on;
    indx_grp={[1:ntime], 1, idx_pad{:}}; indx_grp2=indx_grp;
	for jg=1:ng_ok,
        indx_grp{2}=jg;
        ii=find(~isnan(GRP.data(indx_grp{:})));
        indx_grp2{2}=jg; indx_grp2{1}=ii;
        plot_cfill(GRP.data(indx_grp2{:}),GRP.scales{1}(ii),'o',cols(mod(jg-1,ncol)+1),5); zoom on;
	end
    if nloop>1, pause; end
end
GRP.dims(2)=ng_unique0;
GRP.scales{2}=[1:ng_unique0];
idx_compress={[1:ntime],[1:ng_unique0],indx_dims_full{:}};
GRP.data=GRP.data(idx_compress{:});
GRP.clim=CFAR.scales{jdim2}([1,end]);
%
PARMS.dims(1)=ng_max0;
PARMS.scales{1}=[1:ng_max0];
idx_compress={[1:ng_max0],[1:6],indx_dims_full{:}};
PARMS.data=PARMS.data(idx_compress{:});
PARMS.clim=[min(PARMS.data(:)),max(PARMS.data(:))];
if 0,
    V.ii_img=-1; V.iiq=[1:6]; 
    view5(PARMS,[2 1 3]);
end
