%-------------------------------------------------------------------------% % 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