%-------------------------------------------------------------------------%
%         This software is licensed by Creative Commons BY-NC-SA:         %
%      http://creativecommons.org/licenses/by-nc-sa/3.0/it/legalcode      %
%-------------------------------------------------------------------------%
%
% File: BlockBased3D_Structure(dsites,q,puradius,min_dsites)
%
% Goal: find the data sites located in each of the q^3 cells and in the
%       neighbouring cells
%
% Inputs: dsites:     NX3 matrix representing a set of N data sites
%         q:          number of cells in one direction
%         puradius:   radius of PU subdomains
%         min_dsites: minimum of data sites among two directions
%
% Calls on: countingsort: by B. Moore, from MATLAB Central File Exchange
%
% Outputs: idx_dsites: multiarray containing the indices of the data points 
%                      located in k-th block and in the neighbouring blocks
%
%-------------------------------------------------------------------------%
function [idx_dsites] = BlockBased3D_Structure(dsites,q,puradius,...
    min_dsites)
N = size(dsites,1); k = 1; i1 = 1; % Initialize
% Sort with respect to the first coordinate the data sites
[dsites_sortx,IX] = sortrows(dsites,1);
for k1 = 1:q
    % Build the k1-th strip parallel to the y-axis
    x(k1) = min_dsites  +  puradius*k1; i2 = 1; idx_x = []; % Initialize
    while (i1 <= N) && (dsites_sortx(i1,1) <=x (k1))
        % Find the points located in the k1-th strip parallel to the y-axis
        idx_x(i2) = IX(i1); i1 = i1 + 1; i2 = i2 + 1;
    end
    idx_sortx = countingsort(idx_x,N); % Sort the indices
    % Sort with respect to the second coordinate the data sites
    % located in the k1-th strip parallel to the y-axis
    [dsites_k1x,IX_k1x] = sortrows(dsites(idx_sortx(:),:),2);
    j1 = 1; %Initialize
    for k2 = 1:q
        % Build the k2-th strip parallel to the x-axis
        y(k2) = min_dsites + puradius*k2;
        j2 = 1; idx_y = [];  % Initialize
        while (j1 < i2) && (dsites_k1x(j1,2) <= y(k2))
            % Find the points located in the k2-th strip parallel to the
            % x-axis
            idx_y(j2) = idx_sortx(IX_k1x(j1));
            j1 = j1 + 1; j2 = j2 + 1;
        end
        idx_sorty = countingsort(idx_y,N); % Sort the indices
        % Sort with respect to the third coordinate the data sites
        % located in the k2-th strip parallel to the x-axis
        [dsites_k2y,IX_k2y] = sortrows(dsites(idx_sorty(:),:),3);
        l1 = 1; %Initialize
        for k3 = 1:q
            % Build the k3-th strip parallel to the z-axis
            z(k3) = min_dsites + k3*puradius;
            l2 = 1; idx_z = []; % Initialize
            while (l1 < j2) && (dsites_k2y(l1,3) <=z (k3))
                % Find the points located in the k3-th strip parallel to
                % the z-axis
                idx_z(l2) = idx_sorty(IX_k2y(l1));
                l1 = l1 + 1; l2 = l2 + 1;
            end
            idx_dsites_k{k}  = idx_z';
            k = k + 1;
        end
    end
end
k = 1;
% Find data sites located in the k-th cell and in the neighbouring cells
while k <= q^3
    idx_dsites{k} = idx_dsites_k{k}; % Initialize
    % Reduce the number of neighbouring cells for border cells
    border_cell = 0;
    if k == 1 % Bottom front left corner cell                             
        border_cell = 1; % Reduce the number of neighbouring cells
        idx_dsites{k} = [idx_dsites_k{k};idx_dsites_k{k+1};...
            idx_dsites_k{k+q^2};idx_dsites_k{k+q^2+1};...
            idx_dsites_k{k+q};idx_dsites_k{k+q+1};...
            idx_dsites_k{k+q^2+q};idx_dsites_k{k+q^2+q+1}];
    end
    if k == q^2*(q-1)+1 % Bottom front right corner cell                 
        border_cell = 1; % Reduce the number of neighbouring cells
        idx_dsites{k} = [idx_dsites_k{k};idx_dsites_k{k+1};...
            idx_dsites_k{k-q^2};idx_dsites_k{k-q^2+1};...
            idx_dsites_k{k+q};idx_dsites_k{k+q+1};...
            idx_dsites_k{k-q^2+q};idx_dsites_k{k-q^2+q+1}];
    end
    if k == q % Bottom back left corner cell                            
        border_cell = 1; % Reduce the number of neighbouring cells
        idx_dsites{k} = [idx_dsites_k{k};idx_dsites_k{k-1};...
            idx_dsites_k{k+q^2};idx_dsites_k{k+q^2-1};...
            idx_dsites_k{k+q};idx_dsites_k{k+q-1};...
            idx_dsites_k{k+q^2+q};idx_dsites_k{k+q^2+q-1}];
    end
    if k == (q-1)*q^2+q % Bottom back right corner cell                 
        border_cell = 1; % Reduce the number of neighbouring cells
        idx_dsites{k} = [idx_dsites_k{k};idx_dsites_k{k-1};...
            idx_dsites_k{k-q^2};idx_dsites_k{k-q^2-1};...
            idx_dsites_k{k+q};idx_dsites_k{k+q-1};...
            idx_dsites_k{k-q^2+q};idx_dsites_k{k-q^2+q-1}];
    end
    if k == (q-1)*q+1 % Top front left corner cell                        
        border_cell = 1; % Reduce the number of neighbouring cells
        idx_dsites{k} = [idx_dsites_k{k};idx_dsites_k{k+1};...
            idx_dsites_k{k+q^2};idx_dsites_k{k+q^2+1};...
            idx_dsites_k{k-q};idx_dsites_k{k-q+1};...
            idx_dsites_k{k+q^2-q};idx_dsites_k{k+q^2-q+1}];
    end
    if k == q^3-(q-1) % Top front right corner cell                      
        border_cell = 1; % Reduce the number of neighbouring cells
        idx_dsites{k} = [idx_dsites_k{k};idx_dsites_k{k+1};...
            idx_dsites_k{k-q^2};idx_dsites_k{k-q^2+1};...
            idx_dsites_k{k-q};idx_dsites_k{k-q+1};...
            idx_dsites_k{k-q^2-q};idx_dsites_k{k-q^2-q+1}];
    end
    if k == q^2 % Top back left corner cell                              
        border_cell = 1; % Reduce the number of neighbouring cells
        idx_dsites{k} = [idx_dsites_k{k};idx_dsites_k{k-1};...
            idx_dsites_k{k+q^2};idx_dsites_k{k+q^2-1};...
            idx_dsites_k{k-q};idx_dsites_k{k-q-1};...
            idx_dsites_k{k+q^2-q};idx_dsites_k{k+q^2-q-1}];
    end
    if k == q^3 % Top back right corner cell                             
        border_cell = 1; % Reduce the number of neighbouring cells
        idx_dsites{k} = [idx_dsites_k{k};idx_dsites_k{k-1};...
            idx_dsites_k{k-q^2};idx_dsites_k{k-q^2-1};...
            idx_dsites_k{k-q};idx_dsites_k{k-q-1};...
            idx_dsites_k{k-q^2-q};idx_dsites_k{k-q^2-q-1}];
    end
    % Indices of neighbouring cells                         
    neigh = [k-1,k+1,k+q^2,k+q^2-1,k+q^2+1,k-q^2,k-q^2-1,k-q^2+1,k-q,...
        k-q+1,k-q-1,k-q-q^2,k-q+1-q^2,k-q-1-q^2,k-q+q^2,k-q+1+q^2,...
        k-q-1+q^2,k+q,k+q+1,k+q-1,k+q-q^2,k+q+1-q^2,k+q-1-q^2,k+q+q^2,...
        k+q+1+q^2,k+q-1+q^2];
    for j = 2:q-1                           
        if (k == j) % Bottom left border cells 
            border_cell = 1; % Reduce the number of neighbouring cells
            for i = 1:size(neigh,2)
                if (neigh(i) > 0) && (neigh(i) ~= k-q+q^2+1) && ...
                        (neigh(i) ~= k-q+q^2) && (neigh(i) ~= k-q+q^2-1)
                    idx_dsites{k} = [idx_dsites{k};idx_dsites_k{neigh(i)}];
                end
            end
        end
    end
    for j = 1:q-2                          
        if (k == j*q^2+1) % Bottom front border cells    
                border_cell = 1; % Reduce the number of neighbouring cells
            for i = 1:size(neigh,2)
                if (neigh(i) ~= k-1) && (neigh(i) ~= k+q^2-1) && ...
                        (neigh(i) ~= k-q^2-1) && (neigh(i) ~= k+q-1) && ...
                        (neigh(i) ~= k+q-q^2-1) && ...
                        (neigh(i) ~= k+q+q^2-1) && (neigh(i) ~= k-q) && ...
                        (neigh(i) ~= k-q+1) && (neigh(i) ~= k-q-1) && ...
                        (neigh(i) ~= k-q-q^2+1) && ...
                        (neigh(i) ~= k-q-q^2-1) && ...
                        (neigh(i) ~= k-q-q^2) && ...
                        (neigh(i) ~= k-q+q^2) && ...
                        (neigh(i) ~= k-q+q^2+1) && ...
                        (neigh(i) ~= k-q+q^2-1)
                    idx_dsites{k} = [idx_dsites{k};idx_dsites_k{neigh(i)}];
                end
            end
        end
    end
    for j = 1:q-2                       
        if (k == j*q^2+q) % Bottom back border cells 
            border_cell = 1; % Reduce the number of neighbouring cells
            for i = 1:size(neigh,2)
                if (neigh(i) ~= k+1) && (neigh(i) ~= k+q^2+1) && ...
                        (neigh(i) ~= k-q^2+1) && (neigh(i) ~= k+q+1) && ...
                        (neigh(i) ~= k+q-q^2+1) && ...
                        (neigh(i) ~= k+q+q^2+1) && (neigh(i) ~= k-q) && ...
                        (neigh(i) ~= k-q+1) && (neigh(i) ~= k-q-1) && ...
                        (neigh(i) ~= k-q-q^2+1) && ...
                        (neigh(i) ~= k-q-q^2-1) && ...
                        (neigh(i) ~= k-q-q^2) && ...
                        (neigh(i) ~= k-q+q^2) && ...
                        (neigh(i) ~= k-q+q^2+1) && ...
                        (neigh(i) ~= k-q+q^2-1)
                    idx_dsites{k} = [idx_dsites{k};idx_dsites_k{neigh(i)}];
                end
            end
        end
    end
    for j = 2:q-1          
        if (k == q^2*(q-1)+j) % Bottom right border cells
            border_cell = 1; % Reduce the number of neighbouring cells
            for i = 1:size(neigh,2)
                if (neigh(i) <= q^3) && (neigh(i) ~= k-q) && ...
                        (neigh(i) ~= k-q-1) && (neigh(i) ~= k-q+1) && ...
                        (neigh(i) ~= k-q+1-q^2) && ...
                        (neigh(i) ~= k-q-q^2) && ...
                        (neigh(i) ~= k-q-1-q^2) && ...
                        (neigh(i) ~= k+q^2-q) && ...
                        (neigh(i) ~= k+q^2-q+1) && ...
                        (neigh(i) ~= k+q^2-q-1)
                    idx_dsites{k} = [idx_dsites{k};idx_dsites_k{neigh(i)}];
                end
            end
        end
    end  
    for j1 = 1:q-2                                     
        for j2 = 2:q-1
            if (k == j1*q^2+j2) % Bottom face      
                border_cell = 1; % Reduce the number of neighbouring cells
                for i = 1:size(neigh,2)
                    if  (neigh(i) ~= k-q) && (neigh(i) ~= k-q+1) && ...
                            (neigh(i) ~= k-q-1) && ...
                            (neigh(i) ~= k-q+q^2) && ...
                            (neigh(i) ~= k-q+1+q^2) && ...
                            (neigh(i) ~= k-q-1+q^2) && ...
                            (neigh(i) ~= k-q-q^2) && ...
                            (neigh(i) ~= k-q+1-q^2) && ...
                            (neigh(i) ~= k-q-1-q^2)
                        idx_dsites{k} = [idx_dsites{k}; ...
                            idx_dsites_k{neigh(i)}];
                    end
                end
            end
        end
    end
    for j = 1:q-2          
        if (k == j*q+1) % Front left border cells                                   
            border_cell = 1; % Reduce the number of neighbouring cells
            for i = 1:size(neigh,2)
                if (neigh(i) > 0)  && (neigh(i) ~= k-1) && ...
                        (neigh(i) ~= k-q-1) && (neigh(i) ~= k+q-1) && ...
                        (neigh(i) ~= k+q^2-1) && ...
                        (neigh(i) ~= k-q-1+q^2) && ...
                        (neigh(i) ~= k+q-1+q^2)
                    idx_dsites{k} = [idx_dsites{k};idx_dsites_k{neigh(i)}];
                end
            end
        end
    end   
    for j = 2:q-1                                 
        if (k == j*q) % Back left border cells
            border_cell = 1; % Reduce the number of neighbouring cells
            for i = 1:size(neigh,2)
                if (neigh(i) > 0)  && (neigh(i) ~= k+1) && ...
                        (neigh(i) ~= k-q+1) && (neigh(i) ~= k+q+1) && ...
                        (neigh(i) ~= k+q^2+1) && ...
                        (neigh(i) ~= k-q+1+q^2) && ...
                        (neigh(i) ~= k+q+1+q^2) && ...
                        (neigh(i) ~= k+q+1-q^2)
                    idx_dsites{k} = [idx_dsites{k};idx_dsites_k{neigh(i)}];
                end
            end
        end
    end
    for j = 2:q-1                              
        if (k == (q-1)*q+j) % Top left border cells
            border_cell = 1; % Reduce the number of neighbouring cells
            for i = 1:size(neigh,2)
                if (neigh(i) > 0) && (neigh(i) ~= k+q) && ...
                        (neigh(i) ~= k+q+1) && (neigh(i) ~= k+q-1) && ...
                        (neigh(i) ~= k+q+q^2+1)&& ...
                        (neigh(i) ~= k+q+q^2) && ...
                        (neigh(i) ~= k+q+q^2-1) && ...
                        (neigh(i) ~= k+q+1-q^2) && ...
                        (neigh(i) ~= k+q-q^2) && ...
                        (neigh(i) ~= k+q-q^2-1)
                    idx_dsites{k} = [idx_dsites{k};idx_dsites_k{neigh(i)}];
                end
            end
        end
    end
    for j1 = 2:q-1                                        
        for j2 = 1:q-2
            if (k == j2*q+j1) % Left face   
                border_cell = 1; % Reduce the number of neighbouring cells              
                for i = 1:size(neigh,2)
                    if (neigh(i) > 0)
                        idx_dsites{k} = [idx_dsites{k}; ...
                            idx_dsites_k{neigh(i)}];
                    end
                end
            end
        end
    end    
    for j = 1:q-2                                  
        if (k == (q-1)*q+1+j*q^2) % Top front border cells
            border_cell = 1; % Reduce the number of neighbouring cells
            for i = 1:size(neigh,2)
                if (neigh(i) ~= k-q^2+q+1) && ...
                        (neigh(i) ~= k-q^2+q-1) && ...
                        (neigh(i) ~= k-q^2+q) && (neigh(i) ~= k+q) && ...
                        (neigh(i) ~= k+q+1) && (neigh(i) ~= k+q-1) && ...
                        (neigh(i) ~= k+q+q^2+1)&& ...
                        (neigh(i) ~= k+q+q^2) && ...
                        (neigh(i) ~= k+q+q^2-1) && ...
                        (neigh(i) ~= k+q^2-1) && ...
                        (neigh(i) ~= k-q^2-1) && (neigh(i) ~= k-1) && ...
                        (neigh(i) ~= k-q-1) && ...
                        (neigh(i) ~= k+q^2-1-q) && (neigh(i) ~= k-q^2-1-q)
                    idx_dsites{k} = [idx_dsites{k};idx_dsites_k{neigh(i)}];
                end
            end
        end
    end
    for j = 2:q-1        
        if (k == j*q^2) % Top back border cells 
            border_cell = 1; % Reduce the number of neighbouring cells     
            for i = 1:size(neigh,2)
                if (neigh(i) ~= k-q^2+q+1) && ...
                        (neigh(i) ~= k-q^2+q-1) && ...
                        (neigh(i) ~= k-q^2+q) && (neigh(i) ~= k+q) && ...
                        (neigh(i) ~= k+q+1) && (neigh(i) ~= k+q-1) && ...
                        (neigh(i) ~= k+q+q^2+1)&& ...
                        (neigh(i) ~= k+q+q^2) && ...
                        (neigh(i) ~= k+q+q^2-1) && ...                        
                        (neigh(i) ~= k+q^2+1) && ...
                        (neigh(i) ~= k-q^2+1) && ...
                        (neigh(i) ~= k+1) && (neigh(i) ~= k-q+1) && ...
                        (neigh(i) ~= k+q^2+1-q) && ...
                        (neigh(i) ~= k-q^2+1-q)
                    idx_dsites{k} = [idx_dsites{k};idx_dsites_k{neigh(i)}];
                end
            end
        end
    end
    for j = 2:q-1                                
        if (k == q^3-(q-j)) % Top right border cells  
            border_cell = 1; % Reduce the number of neighbouring cells
            for i = 1:size(neigh,2)
                if (neigh(i) <= q^3) && (neigh(i) ~= k-q^2+q+1) && ...
                        (neigh(i) ~= k-q^2+q-1) && (neigh(i) ~= k-q^2+q)
                    idx_dsites{k} = [idx_dsites{k};idx_dsites_k{neigh(i)}];
                end
            end
        end
    end  
    for j1 = 2:q-1                                           
        for j2 = 1:q-2 
            if (k == (q-1)*q+j1+j2*q^2) % Top face                 
                border_cell = 1; % Reduce the number of neighbouring cells
                for i = 1:size(neigh,2)
                    if (neigh(i) ~= k-q^2+q+1) && ...
                            (neigh(i) ~= k-q^2+q-1) && ...
                            (neigh(i) ~= k-q^2+q) && ...
                            (neigh(i) ~= k+q) && ...
                            (neigh(i) ~= k+q+1) && ...
                            (neigh(i) ~= k+q-1) && ...
                            (neigh(i) ~= k+q+q^2+1)&& ...
                            (neigh(i) ~= k+q+q^2) && ...
                            (neigh(i) ~= k+q+q^2-1)
                        idx_dsites{k} = [idx_dsites{k}; ...
                            idx_dsites_k{neigh(i)}];                       
                    end
                end
            end
        end
    end
    for j = 1:q-2                           
        if (k == q^2*(q-1)+1+(q-(q-j))*q) % Front right border cells                            
            border_cell = 1; % Reduce the number of neighbouring cells           
            for i = 1:size(neigh,2)
                if (neigh(i) <= q^3) && (neigh(i) ~= k-1) && ...
                        (neigh(i) ~= k-q^2-1) && ...
                        (neigh(i) ~= k-q-1) && ...
                        (neigh(i) ~= k-q-1-q^2) && ...
                        (neigh(i) ~= k-q-1+q^2) && ...
                        (neigh(i) ~= k+q-1) && ...
                        (neigh(i) ~= k+q-1-q^2)
                    idx_dsites{k} = [idx_dsites{k};idx_dsites_k{neigh(i)}];
                end
            end
        end
    end   
    for j1 = 1:q-2 
        for j2 = 1:q-2
            if (k == j1*q^2+j2*q+1) % Front face                      
                border_cell = 1; % Reduce the number of neighbouring cells
                for i = 1:size(neigh,2)
                    if (neigh(i) ~= k-q^2-1) && (neigh(i) ~= k-1) && ...
                            (neigh(i) ~= k+q^2-1) && ...
                            (neigh(i) ~= k-q-1) && ...
                            (neigh(i) ~= k-q-1-q^2) && ...
                            (neigh(i) ~= k-q-1+q^2) && ...
                            (neigh(i) ~= k+q-q^2-1)&& ...
                            (neigh(i) ~= k+q-1) && (neigh(i) ~= k+q+q^2-1)
                        idx_dsites{k} = [idx_dsites{k}; ...
                            idx_dsites_k{neigh(i)}];
                    end
                end
            end
        end
    end    
    for j = 1:q-2    
        if (k == q^3-j*q) % Back right border cells                                   
            border_cell = 1; % Reduce the number of neighbouring cells
            for i = 1:size(neigh,2)
                if (neigh(i) <= q^3) && (neigh(i) ~= k+1) && ...
                        (neigh(i) ~= k-q^2+1) && ...
                        (neigh(i) ~= k-q+1) && ...
                        (neigh(i) ~= k-q+1-q^2) && ...
                        (neigh(i) ~= k+q+1) && (neigh(i) ~= k+q+1-q^2)
                    idx_dsites{k} = [idx_dsites{k};idx_dsites_k{neigh(i)}];
                end
            end
        end
    end   
    for j1 = 1:q-2 
        for j2 = 1:q-2 
            if (k == q^2*(q-1)+1+(q-(q-j1))*q+j2) % Right face         
                border_cell = 1; % Reduce the number of neighbouring cells
                for i = 1:size(neigh,2)
                    if  (neigh(i) <= q^3)
                        idx_dsites{k} = [idx_dsites{k}; ...
                            idx_dsites_k{neigh(i)}];
                    end
                end
            end
        end
    end   
    for j1 = 1:q-2 
        for j2 = 2:q-1
            if (k == j1*q^2+q*j2) % Back face                           
                border_cell = 1; % Reduce the number of neighbouring cells
                for i = 1:size(neigh,2)
                    if (neigh(i) ~= k+1) && (neigh(i) ~= k-q^2+1) && ...
                            (neigh(i) ~= k+q^2+1) && ...
                            (neigh(i) ~= k-q+1-q^2) && ...
                            (neigh(i) ~= k-q+1) && ...
                            (neigh(i) ~= k-q+1+q^2) && ...
                            (neigh(i) ~= k-q^2+1+q) && ...
                            (neigh(i) ~= k+q+1) && ...
                            (neigh(i) ~= k+q^2+q+1)
                        idx_dsites{k} = [idx_dsites{k}; ...
                            idx_dsites_k{neigh(i)}];
                    end
                end
            end
        end
    end   
    if  (border_cell == 0)
        for i = 1:size(neigh,2)
            idx_dsites{k} = [idx_dsites{k};idx_dsites_k{neigh(i)}];
        end
    end
    k = k+1;
end