function field = PlaceFieldQuant(pos_map_trial, pos_map_trial_norm, varargin)
% place fields quantification based on
% Mizuseki 2012 Hippocampus and Harvey 2012 Nature
% identification of place cells and quantificatin of fields

% default values
[nTrial, nBin, nUnit, nPermute] = size(pos_map_trial);
min_size = 8;
max_size = 80;
thrd_min_max = 0.3;
half_nTrial = round(nTrial/2);
third_nTrial = round(nTrial/3); % the number of transients should outnumber 1/3 of trial number;
ratio_in_to_out = 3;

% parse input parameters
% Parse parameter list
for i = 1:2:length(varargin),
    if ~ischar(varargin{i}),
		error(['Parameter ' num2str(i+2) ' is not a property']);
    end
    switch(lower(varargin{i}))
        case 'ratio_in_to_out',
			ratio_in_to_out = varargin{i+1};
            if ~isscalar(ratio_in_to_out),
				error('Incorrect value for property ''ratio_in_to_out''');
            end
        case 'min_threshold',
            min_threshold = varargin{i+1};
            if ~isscalar(min_threshold),
                error('Incorrect value for property ''min_threshold''');
            end
		otherwise,
			error(['Unknown property ''' num2str(varargin{i}) '''']);
    end
end

placell_id = zeros(nPermute, nUnit);
field_stab = zeros(1, nUnit);
field_width = nan(1, nUnit);
n_fields = nan(1, nUnit);
field_pdist = nan(1, nUnit);
bin_idx = 1:nBin;

for iUnit = 1:nUnit
    
    for iPermute = 1:nPermute

        pos_map_trial_this = pos_map_trial_norm(:,:,iUnit, iPermute);
        
        if iPermute == 1
            pt_1 = mean(pos_map_trial_this(1:half_nTrial,:),1);
            pt_2 = mean(pos_map_trial_this(half_nTrial+1:end,:),1);
            field_stab(iUnit) = corr(pt_1', pt_2', 'rows', 'pairwise');
        end
        
        pos_tc = mean(pos_map_trial_this(:,:), 1);
        thrd_tmp = (max(pos_tc) - min(pos_tc)) * thrd_min_max;
%         in_field = find(pos_tc-min(pos_tc) > thrd_tmp);
        in_field = pos_tc - min(pos_tc) > thrd_tmp;
        [in_field, n_field, field_n] = FindConsecutive(in_field, min_size, max_size);   

        if ~isempty(in_field)
            if n_field > 1                                                       % account for circular nature
                first_field = field_n == 1;
                last_field = field_n == n_field;
                first_field_idx = in_field(first_field);
                last_field_idx = in_field(last_field);
                if first_field_idx(1) == 1 && last_field_idx(end) == nBin
                    field_n(last_field) = 1;
                    n_field = n_field - 1;
                end
            end
            out_field = bin_idx(~ismember(bin_idx, in_field));
            in_field_resp = pos_tc(in_field);
            out_field_resp = pos_tc(out_field);
            if mean(in_field_resp) > ratio_in_to_out * mean(out_field_resp)
                [~, peak_idx] = max(pos_map_trial_this(:,:), [], 2);
                ok = peak_idx > 1;                                               % exclude silent trials
                peak_idx = peak_idx(ok);

                if sum(ismember(peak_idx, in_field)) >= third_nTrial
                    placell_id(iPermute, iUnit) = 1;
                    pd_mn = 0;
                    for i = 1:n_field                                            % for more than one field
                        ok = ismember(peak_idx, in_field(field_n==i));
                        pd = pdist(peak_idx(ok));
                        ok = pd <= nBin/2;                                       % 100 bins
                        pd(~ok) = nBin-pd(~ok);
                        if ~isempty(pd)
                            pd_mn = pd_mn + mean(pd);
                        end
                    end
                    if iPermute == 1
                        field_width(iUnit) = length(in_field)/n_field;
                        n_fields(iUnit) = n_field;
                        field_pdist(iUnit) = pd_mn / n_field;
                    end
                end
            end
        end
    end

end

field.placell_id = placell_id;
field.stab = field_stab;
field.width = field_width;
field.nfields = n_fields;
field.pdist = field_pdist;       % pairwise drift of field max firing positions across trials

end