function [ih, jh] = gridBasedNonUniformHorizonSearch(coor, delta, ratio, is2cell)
% 这是基于网格方法进行非均匀排布物质点邻域搜索的算法
% 输入参数
% coor ： n×d的坐标矩阵，n表示物质点个数，d表示维度大小
% delta：一个截断半径，也叫horizon size
% ratio：一个大于1的比例参数，当格子划分比较小的时候，增大ratio会相对减少计算用时
% 输出参数
% ih：每个物质点的id
% jh：元胞向量，其中每个元素都是一个向量，代表ih对应的邻域物质带点的id集合
% 在规则网格下，delta=3.015dx情况时，cpu为Intel(R) Xeon(R) CPU E5-2680 v4，耗时为
% 2D: 10000物质点0.15s，100w物质点11.043s
% 3D: 12.5w物质点2.618 秒，100w物质点21.2s
if(nargin==2)  % 如果只是输入两个参数，默认ratio为1
    ratio = 1;
end
if(nargin<=3)
    is2cell = 'no';
end
[np,ndim] = size(coor);                           % 得到物质点个数和维度
vmax      = max(coor);                            % 每个维度的最大值
vmin      = min(coor) - delta * 0.01;             % 每个维度的最小值
ndBin     = floor((vmax-vmin)/(delta*ratio)) + 1; % 每个点在每个维度所处格子区间
nBin      = prod(ndBin);                          % 总共的格子区间个数
hashTable = cell(1,nBin);                         % 每个格子装对应的物质点id
eachSub   = floor((coor - vmin)/(delta*ratio)) + 1;       % 每个物质点对应的格子位置
if(ndim>1)
    subCell   = mat2cell(eachSub, np, ones(1,ndim)); % 把每个物质点所属格子位置分维度存储
    eachIdx   = sub2ind(ndBin, subCell{:});          % 把维度位置换成总的格子指向位置
else
    eachIdx   = eachSub;                             % 如果是1维，就没有必要
end
for i = 1:numel(eachIdx)                             % 对每个格子里面存储物质点id
    j = eachIdx(i);
    hashTable{j} = [hashTable{j}, i];
end
hashHorizon = cell(nBin,1);                          % 初始化horizon可能物质点
for i = 1:nBin                                       % 对每个格子进行循环找到每个格子的邻居格子
    if(isempty(hashTable{i}))                        % 格子里面没有物质点就跳过
        continue;
    end
    if(ndim>1)
        dimCell = cell(1, ndim);
        [dimCell{:}] = ind2sub(ndBin,i);
        for j = 1:ndim
            dimCell{j} = max(1,dimCell{j}-1):min(dimCell{j}+1,ndBin(j));
        end
        [dimCell{:}] = meshgrid(dimCell{:});
        neighborBins = sub2ind(ndBin, dimCell{:}); % 把第i个格子的邻居格子找到
    else
        neighborBins = max(1,i-1):min(i+1,nBin);   % 把第i个格子的邻居格子找到，1维会更简单直接前后两个格子
    end
    p = hashTable(neighborBins(:)');               % 每个格子加上邻近格子的物质点编号
    hashHorizon{i} = [p{:}];                       % 存入hashHorizon中，以方便进一步筛选
end

jh = cell(1,np);
if(ndim>1)                                         % 开始筛选，根据截断半径进行
    delta2 = delta^2;                              % 截断半径delta的平方
    for i = 1:np                                   % 对每个物质点进行循环
        j     = eachIdx(i);                        % 物质点所属的格子位置
        k     = hashHorizon{j};                    % 根据格子位置找到所有可能的邻近物质点
        dist2 = sum((coor(k,:) - coor(i,:)).^2,2); % 计算所有可能邻域物质点到中心物质点的距离平方
        % 问题不大，这里产生的耗时主要来源于减法
        %         dcoor = coor(k,:)-coor(i,:);
        %         dist2 = dcoor(:,1).^2;
        %         for s = 2:ndim
        %             dist2 = dist2  + dcoor(:,s).^2;
        %         end
        q     = dist2>0 & dist2<delta2;             % 满足不重合且在delta范围能选中
        jh{i} = k(q);                               % 存入邻域集合
    end
else
    for i = 1:np                                    % 1维就更简单
        j     = eachIdx(i);
        k     = hashHorizon{j};
        dist  = abs(coor(k,:) - coor(i,:));
        q     = dist>0 & dist<delta;
        jh{i} = k(q);
    end
end
ih = 1:np;

if(strcmp(is2cell,'no'))
    num = cellfun(@numel,jh,'uniform',1);
    ph  = [0, cumsum(num)];
    jh  = [jh{:}]';
    ih  = zeros(size(jh));
    for i = 1:np
        ih(ph(i)+1:ph(i+1)) = i;
    end
end
end