% [profile_x, profile_y, profile_z, final_r_value, lattice] =
% pack_discretization(pack_type, lattice_edge, target_porosity, unit_cells,
% plot_flag)
% creates and discretizes ordered packing of equal spheres arranged in
% body-centered, face-centered, or simple cubic types. Unit cell of a given 
% packing is replicated (unit_cells) times along each Cartesian direction 
% and the resulting geometry is mapped onto a discrete uniform cubic 
% lattice with (lattice_edge)^3 voxels.
%
% EXAMPLE: [p_x, p_y, p_z, r, lattice] =  pack_discretization(1, 50, 0.4, 2, 1);
% 
% Mapping is done (RADIUS_STEPS) times for slight step-wise variation of
% sphere radius around its analytical value to obtain output
% _lattice_ porosity closer to (target_porosity). In some cases, radius
% variation may not change porosity at all.
% 
% Lattice is initialized in (lattice) 3D array with void voxels marked as 0
% and solid - inside spheres - as 1.
%
% If (plot_flag) is not equal to 0 the function draws porosity error vs. 
% sphere radius as well as actual spheres used for discretization.
%
% The code checks output lattice for reflection symmetry along three planes
% and reports number of assymetric voxels for each direction. All three
% numbers have to be 0.
% 
% The function may give warning about twice voxel initialization, 
% which is normal for the case of overlapping spheres.


function [profile_x, profile_y, profile_z, final_r_value, lattice] = pack_discretization(...
  pack_type, ... % 1 for BCC, 2 for FCC, 3 for SC
  lattice_edge, ... % edge of output discrezation domain, in pixels
  target_porosity, ... % desired porosity of the output lattice
  unit_cells, ... % number of unit cell replication along each direction
  plot_flag) % do/not plot porosity errors and packing used for discretization

RADIUS_STEPS = 20; % steps to vary radius in order to better match (target_porosity)
RADIUS_DELTA_PLUS_MINUS = 0.02; % +- 2 percent of voxel size for radius variation

lattice_edge_cube = [lattice_edge, lattice_edge, lattice_edge];
disp('Initialization of...')
switch pack_type
  case 1 % bcc, 8 spheres in case of no periodics
    disp('  BCC packing')
    N = 9;
    packing_sphere = zeros(N, 3);
    
    % eight corner spheres
    for i = 1 : 8
      packing_sphere(i, :) = binary2vector(i - 1, 3) * lattice_edge/unit_cells;
    end
    
    % center sphere
    packing_sphere(N, :) = [lattice_edge/unit_cells/2, lattice_edge/unit_cells/2, lattice_edge/unit_cells/2];
    
    touching_radius = sqrt(3) * lattice_edge/unit_cells / 4; % for the porosity of bcc with touching spheres
  case 2 % fcc, 14 spheres in the case of no periodics
    disp('  FCC packing')
    N = 14;
    packing_sphere = zeros(N, 3);
    
    % eight corner spheres
    for i = 1 : 8
      packing_sphere(i, :) = binary2vector(i - 1, 3) * lattice_edge/unit_cells;
    end

    % three face spheres with one 0 coordinate
    packing_sphere(9:11, :) = unique(perms([0.5 0.5 0]), 'rows') * lattice_edge/unit_cells;
    
    % three face spheres with one lattice_edge coordinate
    packing_sphere(12:14, :) = unique(perms([0.5 0.5 1]), 'rows') * lattice_edge/unit_cells;
    
    touching_radius = sqrt(2) * lattice_edge/unit_cells / 4; % for the porosity of fcc with touching spheres
  case 3 % sc, 1 sphere in the case of no periodics
    disp('  SC packing')
    N = 1;
    packing_sphere = [0.5 0.5 0.5] * lattice_edge/unit_cells;
    
    touching_radius = 0.5 * lattice_edge/unit_cells; % for the porosity of sc with touching spheres
end

indexes_array = zeros(3, lattice_edge_cube(1), lattice_edge_cube(2), lattice_edge_cube(3), 'uint16');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% begin of 'block1'
% this code was extracted from meshgrid routine and modified to reduce 
% memory usage according to specific needs of this program

y = uint16(1:lattice_edge_cube(1));
x = uint16(1:lattice_edge_cube(2));
z = uint16(1:lattice_edge_cube(3));
nx = numel(x);
ny = numel(y);
nz = numel(z);

xx = reshape(x(:),[1 nx 1]);
indexes_array(2, :, :, :) = xx(ones(ny,1),:,ones(nz,1));
clear x xx

yy = reshape(y(:),[ny 1 1]);
indexes_array(1, :, :, :) = yy(:,ones(1,nx),ones(nz,1));
clear y yy

zz = reshape(z(:),[1 1 nz]); % Make sure z is a full page vector.
indexes_array(3, :, :, :) = zz(ones(ny,1),ones(nx,1),:);
clear z zz
% end of 'block1'
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% allocate array for the spheres from all unit cells
spheres = zeros(N*unit_cells^3, 3);

% replicate cells into cube with edge length of 'unit_cells'
for i = 1:unit_cells
  for j = 1:unit_cells
    for k = 1:unit_cells
      start_sphere = (sub2ind([unit_cells, unit_cells, unit_cells], i, j, k) - 1) * N + 1;
      stop_sphere = start_sphere + N - 1;
      spheres(start_sphere:stop_sphere, :) = packing_sphere + repmat([i - 1, j - 1, k - 1]*lattice_edge/unit_cells, N, 1);
    end
  end
end

% remove duplicate spheres and calculate the number of the rest spheres 
spheres = round(spheres*1e6)/1e6;
packing_sphere = unique(spheres, 'rows');
[N, ~] = size(packing_sphere);
disp(['  Total spheres: ', num2str(N)]);

disp('Done')

% modify radii of spheres to obtain target porosity
radius_analyt_scaling_factor = calc_rad_scal_factor(pack_type, target_porosity);
radius_analyt = radius_analyt_scaling_factor * touching_radius;

% define limits and step for the radius variation, which will be done to match 
% the discrete and analytical porosities
r_start = radius_analyt + RADIUS_DELTA_PLUS_MINUS;
r_stop = radius_analyt - RADIUS_DELTA_PLUS_MINUS;

delta_r = (r_start - r_stop)/RADIUS_STEPS;
porosity_errors = zeros(RADIUS_STEPS, 3);

disp([sprintf('\n') 'Looping over ' num2str(RADIUS_STEPS) ' radius values...'])
tic
for k = 1:RADIUS_STEPS
  if mod(k, 2) == 0
   disp([' ' num2str(k) ' out of ' num2str(RADIUS_STEPS)])
  end
  r = r_start - (k-1)*delta_r;
  
  lattice = zeros(lattice_edge, lattice_edge, lattice_edge, 'uint16');
  
  for s = 1 : N % loop over spheres

    indexes_start = max(floor(packing_sphere(s, :) - (r + 1)), 1);
    indexes_stop = min(ceil(packing_sphere(s, :) + (r + 1)), lattice_edge);

    indexes = uint16(squeeze( ((double(indexes_array(1, indexes_start(1):indexes_stop(1), indexes_start(2):indexes_stop(2), indexes_start(3):indexes_stop(3))) - 0.5 - packing_sphere(s, 1)).^2 + ...
                        (double(indexes_array(2, indexes_start(1):indexes_stop(1), indexes_start(2):indexes_stop(2), indexes_start(3):indexes_stop(3))) - 0.5 - packing_sphere(s, 2)).^2 + ...
                        (double(indexes_array(3, indexes_start(1):indexes_stop(1), indexes_start(2):indexes_stop(2), indexes_start(3):indexes_stop(3))) - 0.5 - packing_sphere(s, 3)).^2) < r^2) );

    lattice(indexes_start(1):indexes_stop(1), indexes_start(2):indexes_stop(2), indexes_start(3):indexes_stop(3)) = ...
    lattice(indexes_start(1):indexes_stop(1), indexes_start(2):indexes_stop(2), indexes_start(3):indexes_stop(3)) + indexes;
    clear indexes
  end
  
  total_box_voxels = lattice_edge^3;

  porosity_errors(k,:) = [r, abs((1 - sum(sum(sum(lattice ~= 0)))/total_box_voxels)/target_porosity - 1), 1 - sum(sum(sum(lattice ~= 0)))/total_box_voxels];
  clear dis_array
end
disp('Done')
toc

% if porosity errors are identical, the output radius will be equal
% to its analytical value, otherwise it will be the closest to it
if max(abs(porosity_errors(1:end-1,2)./porosity_errors(2:end,2) - 1)) < 1e-10
  final_r_value = radius_analyt;
  final_r_value_error = porosity_errors(1,2);
else
  % min return first occurence, and radius in 'porosities' decreases
  [por_min, min_location] = min(porosity_errors(:,2));
  if (porosity_errors(min_location, 2) == porosity_errors(min_location + 1, 2) && porosity_errors(min_location, 1) > radius_analyt)
    while (1) 
      min_location = min_location + 1;
      if ~(porosity_errors(min_location + 1, 1) > radius_analyt && porosity_errors(min_location + 1, 2) == por_min)
        break
      end
    end
  end
  
  % if (radius_analyt) is between descrete radii of (equal) minimal porosity 
  % errors, the output radius will be asigned to radius_analyt; otherwise the 
  % output will be the discrete radius value closest to radius_analyt
  if (porosity_errors(min_location + 1, 2) == por_min && radius_analyt < porosity_errors(min_location, 1))
    final_r_value = radius_analyt;
  else
    final_r_value = porosity_errors(min_location, 1);
  end
  final_r_value_error = mean(porosity_errors(min_location, 2));
end

if plot_flag
  figure; hold on; 
  plot(porosity_errors(:, 1), porosity_errors(:, 2), '.-', ...
       [radius_analyt radius_analyt], [0 max(porosity_errors(:, 2))], 'r--', ...
       final_r_value, final_r_value_error, 'go');
  legend('discrete radii', 'analytical radius', 'radius for discretization');
  legend('Location', 'Best');
  xlabel('radius in lattice units');
  ylabel('porosity error');

  figure; hold on; 
  daspect([1 1 1]);
  [x, y, z] = sphere;
  for s = 1 : N
    surf((x*radius_analyt + packing_sphere(s,1)), (y*radius_analyt + packing_sphere(s,2)), (z*radius_analyt + packing_sphere(s,3)))
  end
end

disp([sprintf('\n') 'Discrete (analytical) radius in lattice units: ' num2str(final_r_value) ...
      ' (' num2str(radius_analyt) ')']);

% discretize again, but now with optimal value of r
r = final_r_value;
lattice = zeros(lattice_edge, lattice_edge, lattice_edge, 'uint16');
for s = 1 : N % loop over spheres
  indexes_start = max(floor(packing_sphere(s, :) - (r + 1)), 1);
  indexes_stop = min(ceil(packing_sphere(s, :) + (r + 1)), lattice_edge);

  indexes = squeeze( ((double(indexes_array(1, indexes_start(1):indexes_stop(1), indexes_start(2):indexes_stop(2), indexes_start(3):indexes_stop(3))) - 0.5 - packing_sphere(s, 1)).^2 + ...
                      (double(indexes_array(2, indexes_start(1):indexes_stop(1), indexes_start(2):indexes_stop(2), indexes_start(3):indexes_stop(3))) - 0.5 - packing_sphere(s, 2)).^2 + ...
                      (double(indexes_array(3, indexes_start(1):indexes_stop(1), indexes_start(2):indexes_stop(2), indexes_start(3):indexes_stop(3))) - 0.5 - packing_sphere(s, 3)).^2) < r^2);

  lattice(indexes_start(1):indexes_stop(1), indexes_start(2):indexes_stop(2), indexes_start(3):indexes_stop(3)) = ...
  lattice(indexes_start(1):indexes_stop(1), indexes_start(2):indexes_stop(2), indexes_start(3):indexes_stop(3)) + uint16( indexes );
end

if sum(sum(sum(lattice > 1))) ~= 0
  warning(['  some voxels where initiated twice!'])
end

% solid voxels initiated more than twice will have their values > 1, which
% will affect porosity calulation as well as binary (0 or 1) porosity
% field; let's restore them to unitary value
lattice(lattice > 1) = 1;

porosity_discrete = 1 - sum(sum(sum(lattice)))/total_box_voxels;
disp(['Discrete porosity error relative to its analytical value: ', ...
      num2str(porosity_discrete/target_porosity - 1)]);
profile_x = squeeze( sum(lattice, 1) );
profile_y = squeeze( sum(lattice, 2) );
profile_z = squeeze( sum(lattice, 3) );

% check for symmetry
disp([sprintf('\n') 'Symmetry check (amount of assymetric voxels for each direction):'])
disp([sum(sum(sum( abs(lattice(1:floor(lattice_edge_cube(1)/2), :, :) - lattice(end:-1:ceil(lattice_edge_cube(1)/2+1), :, :) ) ))), ...
      sum(sum(sum( abs(lattice(:, 1:floor(lattice_edge_cube(2)/2), :) - lattice(:, end:-1:ceil(lattice_edge_cube(2)/2+1), :) ) ))), ...
      sum(sum(sum( abs(lattice(:, :, 1:floor(lattice_edge_cube(3)/2)) - lattice(:, :, end:-1:ceil(lattice_edge_cube(3)/2+1)) ) )))] );