%bv_mindistance.m
%For the Harvard Lichtman / Google H01 human cortex EM project
%by Daniel Berger, June 10 2021, adjusted January 19, 2024

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Running this script requires additional files which are provided in the downloadable ZIP package either here:
% https://storage.googleapis.com/h01_paper_public_files/H01_Matlab_analysis_scripts.zip
% or using the link on the H01 landing page.
% To run this script, please download the package, unzip it to your local drive and run the script from its subfolder.
%
% This script computes distances of cell bodies from blood vessels, using the included low-resolution image stacks.
% Saves two Matlab files, bvdist.mat and minbvdist.mat
% minbvdist.mat is used by H01_humananalysis_master.m for the cell type specific blood vessel distance analysis.

bvsourcefilenametemplate='./blood_vessels/goog14_mip4_db_bloodvessels_092_z64_e8.png_s%04d.png';
cbsourcefilenametemplate='./cell_bodies/goog14_cb_400_mip6_export_s%04d.png';

bvsourceslices=0:64:5248;
cbsourceslices=0:128:5248;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The following section defines 'flags' which select which part(s) of this script should be executed.
% Set the flags for the parts you want to execute to 1 and run the script to execute them.

loaddata                   = 1; % Loads the low-resolution image stacks of cell bodies and blood vessels from included files
removenonsurfacepoints     = 1; % Image stack preprocessing
computebvdistancetransform = 1; % Computes a volumetric minimal-distance field of blood vessels and stores to file. Very slow algorithm, can take many hours
findminbvdistance          = 1; % Samples from the volumetric blood vessel distance field at cell body locations and saves to .mat file.
                                % Plotting the distance histograms is done by H01_humananalysis_master.m, section bvmindistanceanalysis

xstretchfactor=11.1/8; % Compensation for section compression caused by ultramicrotome sectioning

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if (loaddata==1)
  for i=1:length(bvsourceslices)
    filename=sprintf(bvsourcefilenametemplate,bvsourceslices(i));
    disp(['Loading ' filename '...']);
    img=imread(filename);
    if (i==1)
      bvv=zeros(size(img,1),size(img,2),length(bvsourceslices),'uint16');
    end
    bvv(:,:,i)=img;
  end

  for i=1:length(cbsourceslices)
    filename=sprintf(cbsourcefilenametemplate,cbsourceslices(i));
    disp(['Loading ' filename '...']);
    img=imread(filename);
    if (i==1)
      cbv=zeros(size(img,1),size(img,2),length(cbsourceslices),'uint16');
    end
    cbv(:,:,i)=img;
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if (removenonsurfacepoints==1)
  disp('Removing interior voxels..');
  a=bvv(2:end-1,2:end-1,2:end-1);
  b=bvv(1:end-2,2:end-1,2:end-1);
  diff=(a~=b);
  b=bvv(3:end,2:end-1,2:end-1);
  diff=bitor(diff,(a~=b));
  b=bvv(2:end-1,1:end-2,2:end-1);
  diff=bitor(diff,(a~=b));
  b=bvv(2:end-1,3:end,2:end-1);
  diff=bitor(diff,(a~=b));
  b=bvv(2:end-1,2:end-1,1:end-2);
  diff=bitor(diff,(a~=b));
  b=bvv(2:end-1,2:end-1,3:end);
  diff=bitor(diff,(a~=b));
  
  a(diff==0)=0;
  bvv(2:end-1,2:end-1,2:end-1)=a;

  a=cbv(2:end-1,2:end-1,2:end-1);
  b=cbv(1:end-2,2:end-1,2:end-1);
  diff=(a~=b);
  b=cbv(3:end,2:end-1,2:end-1);
  diff=bitor(diff,(a~=b));
  b=cbv(2:end-1,1:end-2,2:end-1);
  diff=bitor(diff,(a~=b));
  b=cbv(2:end-1,3:end,2:end-1);
  diff=bitor(diff,(a~=b));
  b=cbv(2:end-1,2:end-1,1:end-2);
  diff=bitor(diff,(a~=b));
  b=cbv(2:end-1,2:end-1,3:end);
  diff=bitor(diff,(a~=b));
  
  a(diff==0)=0;
  cbv(2:end-1,2:end-1,2:end-1)=a;
end

if (computebvdistancetransform==1)
  disp('Computing distance transform, this can take a while (>12 hours)..');
  bvpixelsize=[512 512*xstretchfactor 64*33]; %YXZ of course. stupid Matlab
  d=bwdistsc(bvv,bvpixelsize);
  save('bvdist.mat','d','-v7.3');
end

if (findminbvdistance==1)
  if (~exist('d','var'))
    disp('Loading distance transform..');
    load('bvdist.mat');
  end
  
  disp('Computing distances ...');
  nrcells=max(cbv(:));
  mindist=zeros(nrcells,1);
  d2=d(:,:,1:2:end);
  for i=1:nrcells
    disp(i);
    md=min(d2(cbv==i));
    if (min(size(md))>0)
      mindist(i)=md;
    else
      disp(sprintf('Warning: Cell %d not found.',i));
    end
  end
  save('minbvdist.mat','mindist','-v7.3');
end