%volume_occupancy.m
%For the Harvard Lichtman / Google H01 human cortex EM project
%By Daniel Berger; May 26 2021, January 27, 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 the volume occupancy of the H01 data, in total and per layer (which type of object covers how much of the dataset)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% This section reads the source data into Matlab via VAST and saves it to .mat files.
%%%% Running this requires access to Google Brainmaps and is provided here only for reference.
%%%% You can run the analysis instead using the provided .mat data files.

flags.connecttovast        = 0; % Opens Matlab <-> VAST connection via TCP/IP. Please make sure VAST is running and vasttools.m is running in the same Matlab and can connect
flags.loadlayers           = 0; % Instructs VAST to load the necessary data layers for data extraction
  params.basefolder='D:\Human\H01_Matlab_analysis_scripts\Volume_Occupancy\'; %Please adjust this to where this script and the .vsvr files are on your computer
flags.updategseg           = 0; % Reads the 6-class volume data from Google Brainmaps into Matlab (requires access)
flags.updatetissuemask     = 0; % Reads the approximate tissue mask into Matlab
flags.updatelayermask      = 0; % Reads the cortical layer mask from Google Brainmaps (requires access)
flags.updatelayermasklocal = 0; % Reads the cortical layer mask from VAST segmentation file
flags.updatebvmask         = 0; % Reads the manual blood vessel labels
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Execution control of the different parts of the analysis.
%%%% Set these flags to 1 to execute the respective section of the script.
flags.loaddata         = 1; % Loads the required data into Matlab from included .mat files
flags.analyze          = 1; % Volume occupancy analysis without blood vessels
flags.analyze_with_bvs = 1; % Volume occupancy analysis with blood vessels

params.miplevel=8;
params.imagedatafilename='imagedata_6class.mat';
params.maskdatafilename='maskdata_6class.mat';
params.laydatafilename='laydata_6class.mat';
params.bvdatafilename='bvdata_6class.mat';

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

if (flags.connecttovast==1)
  global vdata;
  if (max(size(vdata))==0)
    disp('WARNING: Not connected to VAST. Please first run VAST, enable the API, then VastTools.m and connect.');
    return;
  end
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  if (vdata.state.isconnected)
    vdata.ui.menu.connectvast.MenuSelectedFcn(); %Disconnect
    vdata.ui.menu.connectvast.MenuSelectedFcn(); %Connect
  else
    vdata.ui.menu.connectvast.MenuSelectedFcn(); %Connect
  end
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  
  vdata.vast.seterrorpopupsenabled(74,0); %Disable error popups when Google sends errors instead of data. VAST will retry
  vdata.vast.seterrorpopupsenabled(75,0); %Disable error popups when Google sends errors instead of data. VAST will retry
end
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if (flags.loadlayers==1)
  disp('Opening layers in VAST...');

  [gseglay,res]=vdata.vast.loadlayer([params.basefolder 'human_SEG_h01_goog14_8nm_6class.vsvr']);
  [llay,res]=vdata.vast.loadlayer([params.basefolder 'human_SEG_h01_goog14_cortical_layers_circular.vsvr']);
  [mslay,res]=vdata.vast.loadlayer([params.basefolder 'goog14_64sec_approximate_6-class_mask.vsseg']);

  info=vdata.vast.getinfo();
  mipscalefactors=vdata.vast.getmipmapscalefactors(0);
end

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

if (flags.updategseg==1)
  miplevel=params.miplevel;
  mipsize=floor(double([info.datasizex info.datasizey info.datasizez])./mipscalefactors(miplevel,:));
  
  disp('Loading data of Google agglomerated segmentation layer...');
  vdata.vast.setselectedlayernr(gseglay);
  gsegimg=zeros(mipsize(2),mipsize(1),mipsize(3),'uint64');
  for i=1:mipsize(3)
    disp(sprintf('Section %d of %d ...',i,mipsize(3)));
    [img,res] = vdata.vast.getemimage(gseglay, params.miplevel, 0, mipsize(1)-1, 0, mipsize(2)-1, i, i);
    gsegimg(:,:,i)=img;
  end
  
  disp(['Saving ' params.imagedatafilename '..']);
  save(params.imagedatafilename,'gsegimg','info','-v7.3');
  disp('  Done');
end

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

if (flags.updatetissuemask==1)
  disp('Loading data of mask segmentation layer...');
  miplevel=params.miplevel;
  mipsize=floor(double([info.datasizex info.datasizey info.datasizez])./mipscalefactors(miplevel,:));
  maskimg=zeros(mipsize(2),mipsize(1),mipsize(3),'uint8');
  for i=1:mipsize(3)
    disp(sprintf('Section %d of %d ...',i,mipsize(3)));
    [segimage, res] = vdata.vast.getsegimageRLEdecoded(params.miplevel, 0, mipsize(1)-1, 0, mipsize(2)-1, i, i, 0,1);
    maskimg(:,:,i)=uint8(segimage);
  end
  
  disp(['Saving ' params.maskdatafilename '..']);
  save(params.maskdatafilename,'maskimg','-v7.3');
  disp('  Done');
end

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

if (flags.updatelayermask==1)
  disp('Loading data of layer layer...');
  miplevel=params.miplevel;
  mipsize=floor(double([info.datasizex info.datasizey info.datasizez])./mipscalefactors(miplevel,:));
 
  vdata.vast.setselectedlayernr(llay);
  layimg=zeros(mipsize(2),mipsize(1),mipsize(3),'uint8');
  for i=1:mipsize(3)
    disp(sprintf('Section %d of %d ...',i,mipsize(3)));
    %simg=vdata.vast.getsegimageRLEdecoded(miplevel, 0, m8size(1)-1, 0, m8size(2)-1, m8sections(i), m8sections(i), 0,1);
    [img,res] = vdata.vast.getemimage(llay, params.miplevel, 0, mipsize(1)-1, 0, mipsize(2)-1, i, i);
    layimg(:,:,i)=uint8(img);
  end

  disp(['Saving ' params.laydatafilename '..']);
  save(params.laydatafilename,'layimg','-v7.3');
  disp('  Done');
end

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

if (flags.updatelayermasklocal==1)
  if (~exist('alaylay','var'))
    [alaylay,res]=vdata.vast.loadlayer([params.basefolder 'H01_alexlayers_jan25_2024_mip8.vsseg']);
  end
  disp('Loading data of cortical layers layer...');
  miplevel=params.miplevel;
  mipscalefactors=vdata.vast.getmipmapscalefactors(0);
  mipsize=floor(double([info.datasizex info.datasizey info.datasizez])./mipscalefactors(miplevel,:));
 
  vdata.vast.setselectedlayernr(alaylay);
  layimg=zeros(mipsize(2),mipsize(1),mipsize(3),'uint8');
  for i=1:mipsize(3)
    disp(sprintf('Section %d of %d ...',i,mipsize(3)));
    [img,res]=vdata.vast.getsegimageRLEdecoded(params.miplevel, 0, mipsize(1)-1, 0, mipsize(2)-1, i, i, 0,1);
    if (res==0) disp('Data loading error.'); end
    layimg(:,:,i)=uint8(img);
  end

  disp(['Saving ' params.laydatafilename '..']);
  save(params.laydatafilename,'layimg','-v7.3');
  disp('  Done');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if (flags.updatebvmask==1)
  if (~exist('bvlay','var'))
    [bvlay,res]=vdata.vast.loadlayer([params.basefolder 'goog14_mip4_db_bloodvessels_092_z64_e20.vsseg']);
  end
  disp('Loading data of blood vessel layer...');
  miplevel=params.miplevel;
  mipscalefactors=vdata.vast.getmipmapscalefactors(0);
  mipsize=floor(double([info.datasizex info.datasizey info.datasizez])./mipscalefactors(miplevel,:));
 
  vdata.vast.setselectedlayernr(bvlay);
  bvimg=zeros(mipsize(2),mipsize(1),mipsize(3),'uint8');
  for i=1:mipsize(3)
    disp(sprintf('Section %d of %d ...',i,mipsize(3)));
    [img,res]=vdata.vast.getsegimageRLEdecoded(params.miplevel, 0, mipsize(1)-1, 0, mipsize(2)-1, i, i, 0,1);
    if (res==0) disp('Data loading error.'); end
    bvimg(:,:,i)=uint8(img>0);
  end

  disp(['Saving ' params.bvdatafilename '..']);
  save(params.bvdatafilename,'bvimg','-v7.3');
  disp('  Done')
end

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

if (flags.loaddata==1)
  disp(['Loading ' params.imagedatafilename '...']);
  load(params.imagedatafilename);
  disp(['Loading ' params.maskdatafilename '...']);
  load(params.maskdatafilename);
  disp(['Loading ' params.laydatafilename '...']);
  load(params.laydatafilename);
  disp(['Loading ' params.bvdatafilename '...']);
  load(params.bvdatafilename);
  
  classnrs=[0 100 101 102 103 104 105 1100 1101 1102 1103];
  gsegimg(gsegimg>1000)=gsegimg(gsegimg>1000)-900;
  gsegimg=uint8(gsegimg);
  
  classnrs=[0 100 101 102 103 104 105 200 201 202 203];
end

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

if (flags.analyze==1)

  labels={'Unlabeled','Axon','Dendrite','Astrocyte','Soma','Cilium','AIS','Myelinated Axon','Myelinated Axon','Fragments 1','Fragments 2'};
  labels2=labels;
  labels3=[];
  totalvol_vox=sum(maskimg(:)>0);
  allpix=gsegimg(maskimg>0);
  
  perc=zeros(11,1);
  for i=1:11
    perc(i)=double(sum(allpix==classnrs(i)))/double(totalvol_vox);
    labels2{i}=[labels{i} sprintf(' (%.2f%%)',perc(i)*100)];
  end
  
  figure(2);
  ex=[1 1 1 1 1 1 1 1 1 1 1];
  pie(perc,ex,labels2);
  
  perc3=zeros(8,1);
  for i=1:7
    perc3(i)=double(sum(allpix==classnrs(i)))/double(totalvol_vox);
    labels3{i}=[labels{i} sprintf(' (%.2f%%)',perc3(i)*100)];
  end
  perc3(8)=double(sum(allpix==classnrs(8))+sum(allpix==classnrs(9))+sum(allpix==classnrs(10))+sum(allpix==classnrs(11)))/double(totalvol_vox);
  labels3{8}=[labels{8} sprintf(' (%.2f%%)',perc3(i)*100)];
  
  figure(3);
  ex=[1 1 1 1 1 1 1 1];
  pie(perc3,ex,labels3);
  
  %%%%%%%%%%%%%%%%%% remove 0
  maskimg2=maskimg;
  maskimg2(gsegimg==0)=0;

  totalvol_vox2=sum(maskimg2(:)>0);
  allpix2=gsegimg(maskimg2>0);
  
  perc=zeros(8,1);
  for i=1:7
    perc(i)=double(sum(allpix2==classnrs(i)))/double(totalvol_vox2);
    labels2{i}=[labels{i} sprintf(' (%.2f%%)',perc(i)*100)];
  end
  perc(8)=double(sum(allpix==classnrs(8))+sum(allpix==classnrs(9))+sum(allpix==classnrs(10))+sum(allpix==classnrs(11)))/double(totalvol_vox2);
  labels2{8}=[labels{8} sprintf(' (%.2f%%)',perc(8)*100)];

  
  figure(4);
  ex=[1 1 1 1 1 1 1];
  piecol=[[0 114 189];[217 83 25];[237 177 32];[126 47 142];[119 172 48];[77 190 238];[162 20 47]]/255;
  ph = pie(perc(2:end),ex,labels2(2:8));
  colormap(piecol);
  %For montage: Width 10, scale font 120. Readjust labels. Manual color transfer from B to A (unknown color mapping)

  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %Analyze by layer
  
  gsegimg2=gsegimg;
  gsegimg2(gsegimg2==100)=1;
  gsegimg2(gsegimg2==101)=2;
  gsegimg2(gsegimg2==102)=3;
  gsegimg2(gsegimg2==103)=4;
  gsegimg2(gsegimg2==104)=5;
  gsegimg2(gsegimg2==105)=6;
  gsegimg2(gsegimg2==200)=7;
  gsegimg2(gsegimg2==201)=8;
  gsegimg2(gsegimg2==202)=9;
  gsegimg2(gsegimg2==203)=10;
  
  laynrvox=zeros(7,1);
  laytypenrvox=zeros(7,10);
  
  for i=1:7
    disp(i);
    maskimg3=maskimg2;
    maskimg3(layimg~=i)=0;
    
    laynrvox(i)=sum(maskimg3(:));
    lvox=gsegimg2.*maskimg3;
    [val,num]=count_unique(lvox(lvox>0));
    laytypenrvox(i,val)=num; 
  end
  
  layperc=laytypenrvox./(laynrvox*ones(1,10));
  figure(5);
  bar(layperc*100,'grouped');
  b=bar(layperc*100,'grouped');
  legend(b,labels(2:end));
  xlabel('Layer');
  ylabel('Volume percent (labeled only)');
  grid on;
  
  layperc2=[layperc(:,1:6) sum(layperc(:,7:10),2)];
  figure(6);
  bar(layperc2*100,'grouped');
  b=bar(layperc2*100,'grouped');
  legend(b,labels(2:end-3));
  ylabel('Volume percent (labeled only)');
  grid on;
  set(gca,'xticklabel',{'Layer 1','Layer 2','Layer 3','Layer 4','Layer 5','Layer 6','White Matter'});
  
  %For montage: Width 10, scale font 120. Readjust labels. Manual color transfer from B to A (unknown color mapping)

  laytypenrvox2=[laytypenrvox(:,1:6) sum(laytypenrvox(:,7:10),2)];
end

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

if (flags.analyze_with_bvs==1)
  
  gsegimg(bvimg>0)=1;
  classnrs=[0 100 101 102 103 104 105 200 201 202 203 1];

  labels={'Unlabeled','Axon','Dendrite','Astrocyte','Soma','Cilium','AIS','Myelinated Axon','Myelinated Axon','Fragments 1','Fragments 2','Blood Vessels'};
  labels2=labels;
  totalvol_vox=sum(maskimg(:)>0);
  allpix=gsegimg(maskimg>0);
  
  perc=zeros(12,1);
  for i=1:12
    perc(i)=double(sum(allpix==classnrs(i)))/double(totalvol_vox);
    labels2{i}=[labels{i} sprintf(' (%.2f%%)',perc(i)*100)];
  end
  
  figure(12);
  ex=[1 1 1 1 1 1 1 1 1 1 1 1];
  pie(perc,ex,labels2);
  
  perc3=zeros(9,1);
  for i=1:7
    perc3(i)=double(sum(allpix==classnrs(i)))/double(totalvol_vox);
    labels3{i}=[labels{i} sprintf(' (%.2f%%)',perc3(i)*100)];
  end
  perc3(8)=double(sum(allpix==classnrs(8))+sum(allpix==classnrs(9))+sum(allpix==classnrs(10))+sum(allpix==classnrs(11)))/double(totalvol_vox);
  labels3{8}=[labels{8} sprintf(' (%.2f%%)',perc3(8)*100)];
  
  perc3(9)=double(sum(allpix==classnrs(12)))/double(totalvol_vox);
  labels3{9}=[labels{12} sprintf(' (%.2f%%)',perc3(9)*100)];
  
  figure(13);
  ex=[1 1 1 1 1 1 1 1 1];
  pie(perc3,ex,labels3);
 

  maskimg2=maskimg;
  maskimg2(gsegimg==0)=0;

  totalvol_vox2=sum(maskimg2(:)>0);
  allpix2=gsegimg(maskimg2>0);
  
  perc=zeros(9,1);
  for i=1:7
    perc(i)=double(sum(allpix2==classnrs(i)))/double(totalvol_vox2);
    labels2{i}=[labels{i} sprintf(' (%.2f%%)',perc(i)*100)];
  end
  perc(8)=double(sum(allpix2==classnrs(8))+sum(allpix2==classnrs(9))+sum(allpix2==classnrs(10))+sum(allpix2==classnrs(11)))/double(totalvol_vox2);
  labels2{8}=[labels{8} sprintf(' (%.2f%%)',perc(8)*100)];

  perc(9)=double(sum(allpix2==classnrs(12)))/double(totalvol_vox2);
  labels2{9}=[labels{12} sprintf(' (%.2f%%)',perc(9)*100)];
  
  figure(14);
  ex=[1 1 1 1 1 1 1 1];
  piecol=[[0 114 189];[217 83 25];[237 177 32];[126 47 142];[119 172 48];[77 190 238];[162 20 47];[255 0 0]]/255;
  ph = pie(perc(2:end),ex,labels2(2:9));
  colormap(piecol);
  %For montage: Width 10, scale font 120. Readjust labels. Manual color transfer from B to A (unknown color mapping)

  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %Analyze by layer
  
  gsegimg2=gsegimg;
  gsegimg2(gsegimg==100)=1;
  gsegimg2(gsegimg==101)=2;
  gsegimg2(gsegimg==102)=3;
  gsegimg2(gsegimg==103)=4;
  gsegimg2(gsegimg==104)=5;
  gsegimg2(gsegimg==105)=6;
  gsegimg2(gsegimg==200)=7;
  gsegimg2(gsegimg==201)=8;
  gsegimg2(gsegimg==202)=9;
  gsegimg2(gsegimg==203)=10;
  gsegimg2(gsegimg==1)=11;
  
  laynrvox=zeros(7,1);
  laytypenrvox=zeros(7,11);
  
  for i=1:7
    disp(i);
    maskimg3=maskimg2;
    maskimg3(layimg~=i)=0;
    
    laynrvox(i)=sum(maskimg3(:));
    lvox=gsegimg2.*maskimg3;
    [val,num]=count_unique(lvox(lvox>0));
    laytypenrvox(i,val)=num;
    
  end
  
  layperc=laytypenrvox./(laynrvox*ones(1,11));
  figure(15);
  bar(layperc*100,'grouped');
  b=bar(layperc*100,'grouped');
  legend(b,labels(2:end));
  xlabel('Layer');
  ylabel('Volume percent (labeled only)');
  grid on;
  
  layperc2=[layperc(:,1:6) sum(layperc(:,7:10),2) layperc(:,11)];
  figure(16);
  bar(layperc2*100,'grouped');
  b=bar(layperc2*100,'grouped');
  legend(b,labels([2:end-4 end]));
  ylabel('Volume percent (labeled only)');
  grid on;
  set(gca,'xticklabel',{'Layer 1','Layer 2','Layer 3','Layer 4','Layer 5','Layer 6','White Matter'});
  
  %For montage: Width 10, scale font 120. Readjust labels. Manual color transfer from B to A (unknown color mapping)

  laytypenrvox2=[laytypenrvox(:,1:6) sum(laytypenrvox(:,7:10),2) laytypenrvox(:,11)];
  
  disp('Volume percent per layer (labeled only):');
  disp(100*layperc2');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The 'count_unique' function is not a core function of Matlab; we are using the code below

function [uniques,numUnique] = count_unique(x,option)
%COUNT_UNIQUE  Determines unique values, and counts occurrences
%   [uniques,numUnique] = count_unique(x)
%
%   This function determines unique values of an array, and also counts the
%   number of instances of those values.
%
%   This uses the MATLAB builtin function accumarray, and is faster than
%   MATLAB's unique function for intermediate to large sizes of arrays for integer values.  
%   Unlike 'unique' it cannot be used to determine if rows are unique or 
%   operate on cell arrays.
%
%   If float values are passed, it uses MATLAB's logic builtin unique function to
%   determine unique values, and then to count instances.
%
%   Descriptions of Input Variables:
%   x:  Input vector or matrix, N-D.  Must be a type acceptable to
%       accumarray, numeric, logical, char, scalar, or cell array of
%       strings.
%   option: Acceptable values currently only 'float'.  If 'float' is
%           specified, the input x vector will be treated as containing
%           decimal values, regardless of whether it is a float array type.
%
%   Descriptions of Output Variables:
%   uniques:    sorted unique values
%   numUnique:  number of instances of each unique value
%
%   Example(s):
%   >> [uniques] = count_unique(largeArray);
%   >> [uniques,numUnique] = count_unique(largeArray);
%
%   See also: unique, accumarray

% Author: Anthony Kendall
% Contact: anthony [dot] kendall [at] gmail [dot] com
% Created: 2009-03-17

testFloat = false;
if nargin == 2 && strcmpi(option,'float')
    testFloat = true;
end

nOut = nargout;
if testFloat
    if nOut < 2
        [uniques] = float_cell_unique(x,nOut);
    else
        [uniques,numUnique] = float_cell_unique(x,nOut);
    end
else
    try %this will fail if the array is float or cell
        if nOut < 2
            [uniques] = int_log_unique(x,nOut);
        else
            [uniques,numUnique] = int_log_unique(x,nOut);
        end
    catch %default to standard approach
        if nOut < 2
            [uniques] = float_cell_unique(x,nOut);
        else
            [uniques,numUnique] = float_cell_unique(x,nOut);
        end
    end
end

end

function [uniques,numUnique] = int_log_unique(x,nOut)
%First, determine the offset for negative values
minVal = min(x(:));

%Check to see if accumarray is appropriate for this function
maxIndex = max(x(:)) - minVal + 1;
if maxIndex / numel(x) > 1000
    error('Accumarray is inefficient for arrays when ind values are >> than the number of elements')
end

%Now, offset to get the index
index = x(:) - minVal + 1;

%Count the occurrences of each index value
numUnique = accumarray(index,1);

%Get the values which occur more than once
uniqueInd = (1:length(numUnique))';
uniques = uniqueInd(numUnique>0) + minVal - 1;

if nOut == 2
    %Trim the numUnique array
    numUnique = numUnique(numUnique>0);
end
end 

function [uniques,numUnique] = float_cell_unique(x,nOut)

if ~iscell(x)
    %First, sort the input vector
    x = sort(x(:));
    numelX = numel(x);
    
    %Check to see if the array type needs to be converted to double
    currClass = class(x);
    isdouble = strcmp(currClass,'double');
    
    if ~isdouble
        x = double(x);
    end
    
    %Check to see if there are any NaNs or Infs, sort returns these either at
    %the beginning or end of an array
    if isnan(x(1)) || isinf(x(1)) || isnan(x(numelX)) || isinf(x(numelX))
        %Check to see if the array contains nans or infs
        xnan = isnan(x);
        xinf = isinf(x);
        testRep = xnan | xinf;
        
        %Remove all of these from the array
        x = x(~testRep);
    end
    
    %Determine break locations of unique values
    uniqueLocs = [true;diff(x) ~= 0];
else
    isdouble = true; %just to avoid conversion on finish
    
    %Sort the rows of the cell array
    x = sort(x(:));
    
    %Determine unique location values
    uniqueLocs = [true;~strcmp(x(1:end-1),x(2:end)) ~= 0] ;
end

%Determine the unique values
uniques = x(uniqueLocs);

if ~isdouble
    x = feval(currClass,x);
end

%Count the number of duplicate values
if nOut == 2
    numUnique = diff([find(uniqueLocs);length(x)+1]);
end
end