Add structure tensor code

BlockIndices.m

+function [B,E,BB,BE] = BlockIndices(VolSiz, BlkSiz, BdrSiz)
+% BLOCKINDICES - Get block indices used for block processing
+%
+% [B,E,BB,EB] = BlockIndices(VolSiz, BlkSiz)
+% [B,E,BB,EB] = BlockIndices(VolSiz, BlkSiz, BdrSiz)
+%
+% INPUT
+% VolSiz - N-D vector with the size of the volume
+% BlkSiz - Scalar or N-D vector with the size of the blocks
+% BdrSiz - (Optional) Scalar or N-D vector with the size of the border
+%
+% OUTPUT
+% B      - N-D Cell with block start points
+% E      - N-D Cell with block end points
+% BB     - N-D Cell with block start points with borders. If BdrSiz == 0,
+%          this is equal to B
+% EB     - N-D Cell with block end points with borders. If BdrSiz == 0, 
+%          this is equal to E
+%
+% Patrick M. Jensen, 2019, ROCKWOOL
+
+% Handle optional
+if nargin < 3, BdrSiz = 0; end
+
+% Process inputs
+N = length(VolSiz);
+
+if length(BlkSiz) == 1
+    BlkSiz = repmat(BlkSiz,1,N);
+elseif length(BlkSiz) ~= N
+    error('BlkSiz must be a scalar, or have same length as VolSiz');
+end
+if length(BdrSiz) == 1
+    BdrSiz = repmat(BdrSiz,1,N);
+elseif length(BlkSiz) ~= N
+    error('If given, BdrSiz must be a scalar, or have same length as VolSiz');
+end
+
+
+B = cell(1,N);
+E = cell(1,N);
+BB = cell(1,N);
+BE = cell(1,N);
+for i = 1:N
+    % Block start points
+    B{i} = 1:BlkSiz(i):VolSiz(i);
+    
+    % Block end points
+    E{i} = min(B{i} + BlkSiz(i) - 1, VolSiz(i));
+    
+    % Block start points with borders
+    BB{i} = max(1, B{i} - BdrSiz(i));
+    
+    % Block end points with borders
+    BE{i} = min(E{i} + BdrSiz(i), VolSiz(i));
+end

CalcEigRealSymm3.m

+function varargout = CalcEigRealSymm3(A11,A22,A33,A12,A13,A23,nval,nvec)
+% CALCEIGREALSYMM3 Eigenvalues and -vectors for symmetric real matrices
+% [...] = CalcEigRealSymm3(A11,A22,A33,A12,A13,A23,NVal,NVec)
+%
+% NOTE: This function is mostly meant for internal use by EIGREALSYMM3 and
+% STRUCTURETENSOREIG3.
+%
+% INPUT
+% A11,A22,A33,A12,A13,A23 - N x M x L arrays with matrix elements
+% NVal                    - Return N smallest eigenvalues
+% NVec                    - Return eigenvectors for N smallest eigenvalues
+%
+% OUTPUT
+% [...] - NVal arrays with eigenvalues followed by 3*NVec arrays with
+%         eigenvector componentes (in order X1,Y1,Z1,X2,Y2,Z2,X3,Y3,Z3).
+%         All arrays have size N x M x L.
+%
+% See also EIGREALSYMM3 STRUCTURETENSOREIG3
+%
+% Patrick M. Jensen, 2019, ROCKWOOL
+
+% Implements the algorithm in:
+%
+% Scherzinger, W. M; Dohrman, C. R
+% A robust algorithm for finding the eigenvalues and eigenvectors of 3 x 3
+% symmetric matrices
+% Comput. Methods Appl. Mech. Engrg., 2008
+% DOI: 10.1016/j.cma.2008.03.031
+
+% NOTE: This code is optimized for speed, *not* readability. To understand
+% each step, please read the above paper first!
+
+% Save original size to reshape at the end
+origsiz = size(A11);
+
+varargout = cell(1,nval+3*nvec);
+
+% Reshape all inputs to 1D arrays
+A11 = reshape(A11,[],1);
+A22 = reshape(A22,[],1);
+A33 = reshape(A33,[],1);
+A12 = reshape(A12,[],1);
+A13 = reshape(A13,[],1);
+A23 = reshape(A23,[],1);
+
+siz = size(A11);
+
+% Convert to deviatoric matrix
+trA = (A11 + A22 + A33)/3;
+A11 = A11 - trA;
+A22 = A22 - trA;
+A33 = A33 - trA;
+
+% Find the most distinct eigenvalue
+Eta1 = MostDistEigVal(A11,A22,A33,A12,A13,A23);
+
+% Find the most distinct eigenvector v1, and s1 and s2
+[V1,S1,S2] = MostDistEigBasis(A11,A22,A33,A12,A13,A23,Eta1);
+
+% Find the other eigenvalues
+[Eta2,Eta3] = OtherEigVals(A11,A22,A33,A12,A13,A23,S1,S2);
+
+% Find the other eigenvectors
+W1 = ComputeW1(A11,A22,A33,A12,A13,A23,S1,S2,Eta2);
+clear S1 S2;
+V2 = RowCross(V1,W1);
+clear W1;
+V2 = bsxfun(@rdivide, V2, sqrt(RowNorms2(V2)));
+V3 = RowCross(V2,V1);
+
+% Get the smallest eigenvalue and corresponding eigencectors
+% NOTE: Since we know there are 3 eigenvalues, this is way of sorting is
+% faster than calling the sort function.
+AllEta = [Eta1,Eta2,Eta3];
+Es = zeros(size(AllEta),'like',AllEta);
+Idx = zeros(size(AllEta),'like',AllEta);
+[Es(:,1),Idx(:,1)] = min(AllEta,[],2);
+[Es(:,3),Idx(:,3)] = max(AllEta,[],2);
+Idx(:,3) = Idx(:,3) + (Es(:,1) == Es(:,3));
+Idx(:,2) = 6-Idx(:,1)-Idx(:,3);
+for i = 1:3
+    Es(Idx(:,2) == i,2) = AllEta(Idx(:,2) == i,i);
+end
+
+clear Eta2 Eta3;
+Es = Es + repmat(trA,1,3);
+for i = 1:nval
+    varargout{i} = reshape(Es(:,i),origsiz);
+end
+
+for i = 1:nvec
+    Xi = zeros(siz,'like',A11);
+    Yi = zeros(siz,'like',A11);
+    Zi = zeros(siz,'like',A11);
+    
+    Xi(Idx(:,i) == 1) = V1(Idx(:,i) == 1,1);    
+    Xi(Idx(:,i) == 2) = V2(Idx(:,i) == 2,1);
+    Xi(Idx(:,i) == 3) = V3(Idx(:,i) == 3,1);
+    Xi(Eta1 == 0) = 1;
+    
+    Yi(Idx(:,i) == 1) = V1(Idx(:,i) == 1,2);
+    Yi(Idx(:,i) == 2) = V2(Idx(:,i) == 2,2);
+    Yi(Idx(:,i) == 3) = V3(Idx(:,i) == 3,2);
+    Yi(Eta1 == 0) = 0;
+    
+    Zi(Idx(:,i) == 1) = V1(Idx(:,i) == 1,3);
+    Zi(Idx(:,i) == 2) = V2(Idx(:,i) == 2,3);
+    Zi(Idx(:,i) == 3) = V3(Idx(:,i) == 3,3);
+    Zi(Eta1 == 0) = 0;
+    
+    varargout{nval+(i-1)*3+1} = reshape(Xi,origsiz);
+    varargout{nval+(i-1)*3+2} = reshape(Yi,origsiz);
+    varargout{nval+(i-1)*3+3} = reshape(Zi,origsiz);
+end
+end
+
+function W1 = ComputeW1(A11,A22,A33,A12,A13,A23,S1,S2,Eta2)
+% Helper function which returns w1
+
+% Save the rows of A - Eta2*I
+AEta1 = [A11 - Eta2, A12, A13];
+AEta2 = [A12, A22 - Eta2, A23];
+AEta3 = [A13, A23, A33 - Eta2];
+
+% Compute u1 and u2, select the largest and normalize
+U1 = [dot(AEta1,S1,2), dot(AEta2,S1,2), dot(AEta3,S1,2)];
+U2 = [dot(AEta1,S2,2), dot(AEta2,S2,2), dot(AEta3,S2,2)];
+clear AEta1 AEta2 AEta3;
+[LU1s, LIdx] = max([RowNorms2(U1), RowNorms2(U2)],[],2);
+
+% We use bsxfun instead of logical indexing because it is faster. The
+% reason is likely that this way we don't indexing (which means branching)
+% and everything is the same operations (which is especially good on GPUs).
+W1 = bsxfun(@times, U1, LIdx == 1);
+W1 = W1 + bsxfun(@times, U2, LIdx == 2);
+clear U1 U2;
+W1 = bsxfun(@rdivide, W1, sqrt(LU1s));
+% In case both u1 and u2 were zero vectors they were both eigenvectors
+% for eta2, so returning s1 will make the rest of the algorithm work
+W1(abs(LU1s) < eps,:) = S1(abs(LU1s) < eps,:);
+end
+
+function [Eta2, Eta3] = OtherEigVals(A11,A22,A33,A12,A13,A23,S1,S2)
+% Helper function which returns the other eigenvalues eta2 and eta3
+
+[Ap22,Ap33,Ap23] = SBasisLowRight(A11,A22,A33,A12,A13,A23,S1,S2);
+D = Ap22 - Ap33;
+trAp = Ap22 + Ap33;
+SgnD = -1 + 2 * (D > 0);
+Eta2 = 0.5*(trAp - SgnD.*sqrt(D.^2 + 4*Ap23.^2));
+Eta3 = trAp - Eta2;
+end
+
+function [Ap11,Ap22,Ap12] = SBasisLowRight(A11,A22,A33,A12,A13,A23,S1,S2)
+% Helper function which returns the lower right 2x2 part of the deviatoric
+% matrix in the (v1,s1,s2) basis
+
+% Save the rows of A
+A1 = [A11,A12,A13];
+A2 = [A12,A22,A23];
+A3 = [A13,A23,A33];
+
+% Compute matrix components
+Ap11 = dot(S1, [dot(A1,S1,2), dot(A2,S1,2), dot(A3,S1,2)],2);
+Ap22 = dot(S2, [dot(A1,S2,2), dot(A2,S2,2), dot(A3,S2,2)],2);
+Ap12 = dot(S1, [dot(A1,S2,2), dot(A2,S2,2), dot(A3,S2,2)],2);
+end
+
+function [V1,S1,S2] = MostDistEigBasis(A11,A22,A33,A12,A13,A23,Eta)
+% Helper function which returns the most distinct eigenvector, s1, and s2
+
+% Compute the basis vector projections
+R1 = [A11 - Eta, A12, A13];
+R2 = [A12, A22 - Eta, A23];
+R3 = [A13, A23, A33 - Eta];
+
+% Pick the largest r vector and make sorted vectors
+[LR1s, RIdx] = max([RowNorms2(R1),RowNorms2(R2),RowNorms2(R3)],[],2);
+
+% We use bsxfun instead of logical indexing because it is faster. The
+% reason is likely that this way we don't indexing (which means branching)
+% and everything is the same operations (which is especially good on GPUs).
+R1s = bsxfun(@times, R1, RIdx == 1);
+R1s = R1s + bsxfun(@times, R2, RIdx == 2);
+R1s = R1s + bsxfun(@times, R3, RIdx == 3);
+
+R2s = bsxfun(@times, R1, RIdx == 2);
+R2s = R2s + bsxfun(@times, R2, RIdx ~= 2);
+
+R3s = bsxfun(@times, R1, RIdx == 3);
+R3s = R3s + bsxfun(@times, R3, RIdx ~= 3);
+
+clear R1 R2 R3 RIdx;
+
+% Compute the s1, t1 and t2 vectors
+S1 = bsxfun(@rdivide, R1s, sqrt(LR1s));
+clear LR1s;
+T2 = R2s - bsxfun(@times, dot(S1,R2s,2), S1);
+T3 = R3s - bsxfun(@times, dot(S1,R3s,2), S1);
+
+% Select the largest of the t1 and t2 vectors, normalize and compute the
+% v1 vectors
+[LT2s, TIdx] = max([RowNorms2(T2),RowNorms2(T3)],[],2);
+
+% We use bsxfun instead of logical indexing because it is faster. The
+% reason is likely that this way we don't indexing (which means branching)
+% and everything is the same operations (which is especially good on GPUs).
+S2 = bsxfun(@times, T2, TIdx == 1);
+S2 = S2 + bsxfun(@times, T3, TIdx == 2);
+S2 = bsxfun(@rdivide, S2, sqrt(LT2s));
+
+V1 = RowCross(S1,S2);
+end
+
+function Eta = MostDistEigVal(A11,A22,A33,A12,A13,A23)
+% Helper function which returns the most distinct eigenvalue
+
+P = sqrt((A11.^2 + A22.^2 + A33.^2 + 2*(A12.^2 + A13.^2 + A23.^2))/6);
+J3Scaled = 0.5*Det3(A11./P,A22./P,A33./P,A12./P,A13./P,A23./P);
+
+% Small numerical errors might cause the argument for acos to have an
+% abslute value greater than 1, which results in a complex result.
+% Therefore, to avoid this, we just truncate it prior to computation.
+Alpha = acos(max(-1, min(1, J3Scaled)))/3;
+clear J3Scaled;
+
+Idx = Alpha > pi/6;
+Eta = 2*P.*cos(Alpha + Idx.*2*pi/3);
+Eta(isnan(Eta)) = 0;
+end
+
+function D = Det3(A11,A22,A33,A12,A13,A23)
+% Helper function which returns the determinant of a 3x3 symetric matrix
+
+D = A11.*A22.*A33 - A11.*A23.^2 - A33.*A12.^2 + 2*A12.*A13.*A23 - ...
+    A13.^2.*A22;
+end
+
+function V = RowCross(V1,V2)
+% Helper function which returns cross product of vectors stacked in rows
+
+% For some reason, it is *much* faster to this manually, even though it is
+% the same as the cross function
+V = [V1(:,2).*V2(:,3) - V1(:,3).*V2(:,2),...
+     V1(:,3).*V2(:,1) - V1(:,1).*V2(:,3),...
+     V1(:,1).*V2(:,2) - V1(:,2).*V2(:,1)];
+end
+
+function N = RowNorms2(V)
+% Helper function which returns the square norms of vectors stacked in rows
+
+N = sum(V.^2,2);
+end

CalcStructureTensor3.m

+function [T11,T22,T33,T12,T13,T23] = CalcStructureTensor3(V,sigma,rho)
+% CALCSTRUCTURETENSOR3 Compute the structure tensors for a 3D volume
+% [T11,T22,T33,T12,T13,T23] = CalcStructureTensor3(V,sigma,rho)
+%
+% NOTE: This function is mostly meant for internal use by STRUCTURETENSOR3
+% and STRUCTURETENSOREIG3.
+%
+% INPUT
+% V     - Input 3D volume of size N x M x L
+% sigma - Std. dev. of convolution kernel for derivatives
+% rho   - Std. dev. of convolution kernel for structure tensor
+%
+% OUTPUT
+% T11, T22, T33, T12 T13, T23 - N x M x L arrays with the elements of the
+%                               structure tensor
+% Patrick M. Jensen, 2019, ROCKWOOL
+
+% Based on the method described in
+%
+% Krause, M.; et. al.
+% Determination of the fibre orientation in composites using the structure
+% tensor and local X-ray transform
+% J. Mater. Sci., 2010
+% DOI: 10.1007/s10853-009-4016-4
+
+% Gaussian derivatives
+x = (-ceil(2*sigma) : ceil(2*sigma))';
+g = 1/(sigma*sqrt(2*pi))*exp(-x.^2/(2*sigma.^2));
+g = g / sum(g);
+dg = -x.*g/sigma.^2;
+Vx = imfilter(imfilter(imfilter(V,dg,'replicate'),g','replicate'),permute(g,[2,3,1]),'replicate');
+Vy = imfilter(imfilter(imfilter(V,dg','replicate'),g,'replicate'),permute(g,[2,3,1]),'replicate');
+Vz = imfilter(imfilter(imfilter(V,permute(dg,[2,3,1]),'replicate'),g,'replicate'),g','replicate');
+
+% Compute structure tensor
+stfsiz = 2*ceil(2*rho)+1;
+SmoothST = @(V) imgaussfilt3(V,rho,'FilterSize',stfsiz,...
+    'FilterDomain','spatial');
+
+T11 = SmoothST(Vx.^2);
+T12 = SmoothST(Vx.*Vy);
+T13 = SmoothST(Vx.*Vz);
+clear Vx;
+T22 = SmoothST(Vy.^2);
+T23 = SmoothST(Vy.*Vz);
+clear Vy;
+T33 = SmoothST(Vz.^2);
+end

EigRealSymm3.m

+function [E,X,Y,Z] = EigRealSymm3(A11,A22,A33,A12,A13,A23,varargin)
+% EIGREALSYMM3 Eigenvalues and -vectors for symmetric real 3 x 3 matrices
+% [E,X,Y,Z] = EigRealSymm3(A11,A22,A33,A12,A13,A23)
+% [E,X,Y,Z] = EigRealSymm3(___,Name,Value,...)
+%
+% INPUT
+% A11, A22, A33, A12, A13, A23 - N x M x L arrays with matrix elements
+%
+% KEY-VALUE PARAMETERS
+% NumEigVals - Return N smallest eigenvalues (default: 3)
+% NumEigVecs - Return eigenvectors for N smallest eigenvalues (default: 3)
+% Verbose    - Print status information (default: false)
+% UseGpu     - Use the GPU. Must be logical or text with 'auto'. If value
+%              is 'auto' the GPU is used if any is available 
+%              (default: 'auto')
+% BlockSize  - Block size for block processing (default: size(V))
+%
+% OUTPUT
+% E - NumEigVals x 1 cell with N x M x L arrays with eigenvalues
+% X - NumEigVecs x 1 cell with N x M x L array with the x-components
+% Y - NumEigVecs x 1 cell with N x M x L array with the y-components
+% Z - NumEigVecs x 1 cell with N x M x L array with the z-components
+%
+% Patrick M. Jensen, 2019, ROCKWOOL
+
+% Handle inputs
+narginchk(6,inf);
+if ~isequal(size(A11), size(A22), size(A33),...
+        size(A12), size(A13), size(A23))
+    error('Input arrays must have the same size');
+end
+validateInputArray = @(x) validateattributes(x,...
+    {'double','single'},{'real','nonsparse'});
+prs = inputParser;
+% Required
+addRequired(prs,'A11',validateInputArray);
+addRequired(prs,'A22',validateInputArray);
+addRequired(prs,'A33',validateInputArray);
+addRequired(prs,'A12',validateInputArray);
+addRequired(prs,'A13',validateInputArray);
+addRequired(prs,'A23',validateInputArray);
+% Key-value parameters
+addParameter(prs,'Verbose',false,@(x) validateattributes(x,...
+    {'logical'},{'scalar'}));
+addParameter(prs,'UseGpu','auto',@(x) validateattributes(x,...
+    {'logical','string','char'},{}));
+addParameter(prs,'BlockSize',size(A11),@(x) validateattributes(x,...
+    {'numeric'},{'integer','vector','positive'}));
+addParameter(prs,'NumEigVals',3,@(x) validateattributes(x,...
+    {'numeric'},{'integer','nonnegative','<=',3,'scalar'}));
+addParameter(prs,'NumEigVecs',3,@(x) validateattributes(x,...
+    {'numeric'},{'integer','nonnegative','<=',3,'scalar'}));
+parse(prs,A11,A22,A33,A12,A13,A23,varargin{:});
+Opts = prs.Results;
+switch validatestring(string(Opts.UseGpu),{'false','true','auto'})
+    case 'false', Opts.UseGpu = false;
+    case 'true', Opts.USeGpu = true;
+    case 'auto', Opts.UseGpu = gpuDeviceCount > 0;
+end
+if Opts.UseGpu && gpuDeviceCount == 0
+    error('No GPU available');
+end
+
+if Opts.Verbose, fprintf('Doing diagonalization\n'); end
+NVal = Opts.NumEigVals;
+NVec = Opts.NumEigVecs;
+% Send off for block processing
+% TODO: Don't use volume blocks, as we don't care about spatial location.
+% This is likely not cache friendly!
+nout = NVal+3*NVec;
+Res = cell(1,nout);
+[Res{:}] = Gpu3dBlockProc({A11,A22,A33,A12,A13,A23},Opts.BlockSize,...
+    @(A11,A22,A33,A12,A13,A23) ...
+    CalcEigRealSymm3(A11,A22,A33,A12,A13,A23,NVal,NVec),...
+    'Verbose',Opts.Verbose,'UseGpu',Opts.UseGpu,'NumOut',nout);
+% Allocate output cells
+E = cell(1,Opts.NumEigVals);
+X = cell(1,Opts.NumEigVecs);
+Y = cell(1,Opts.NumEigVecs);
+Z = cell(1,Opts.NumEigVecs);
+% Gather up results
+for i = 1:NVal
+    E{i} = Res{i};
+end
+for i = 1:NVec
+    X{i} = Res{NVal+(i-1)*3+1};
+    Y{i} = Res{NVal+(i-1)*3+2};
+    Z{i} = Res{NVal+(i-1)*3+3};
+end
+if Opts.Verbose, fprintf('Done\n'); end
+
+end

Gpu3dBlockProc.m

+function varargout = Gpu3dBlockProc(Vi,BlkSz,Fun,varargin)
+% GPU3DBLOCKPROC GPU block processing for 3D volumes
+% Vo = Gpu3dBlockProc(Vi,BlkSiz,Fun)
+% Vo = Gpu3dBlockProc(___,Name,Value,...)
+%
+% INPUT
+% Vi    - Volume or cell with NumIn input volumes
+% BlkSz - Block size - must be [m n l] vector or a scalar
+% Fun   - Function to apply - must take NumIn inputs and return NumOut
+%         outputs
+%
+% KEY-VALUE PARAMETERS
+% BorderSize - [n m l] vector or scalar with the size of the block border 
+%              (default: 0)
+% Verbose    - Print status messages (default: false)
+% NumOut     - Number of outputs for Fun (default: 1)
+% UseGpu     - Do processing on the GPU (default: true)
+%
+% OUTPUT
+% [...] - NumOut x 1 (varargout) cell with output volumes
+%
+% See also BLOCKINDICES
+%
+% Patrick M. Jensen, 2019, ROCKWOOL
+
+% This function was heavily inspired by MATLABs blockproc function, as
+% well as:
+%
+% blockproc3 from the Gerardus toolbox
+% https://github.com/vigente/gerardus/blob/master/matlab/FiltersToolbox/blockproc3.m
+
+% Handle inputs
+narginchk(3,inf);
+prs = inputParser;
+% Required
+addRequired(prs,'Vi',@(x) validateattributes(x,...
+    {'cell'},{'vector','nonempty'}));
+addRequired(prs,'BlkSz',@(x) validateattributes(x,...
+    {'numeric'},{'integer','vector','positive'}));
+addRequired(prs,'Fun',@(x) validateattributes(x,...
+    {'function_handle'},{}));
+% Key-value parameters
+addParameter(prs,'BorderSize',0,@(x) validateattributes(x,...
+    {'numeric'},{'vector','vector','nonnegative'}));
+addParameter(prs,'Verbose',false,@(x) validateattributes(x,...
+    {'logical'},{'scalar'}));
+addParameter(prs,'NumOut',1,@(x) validateattributes(x,...
+    {'numeric'},{'scalar','positive','integer'}));
+addParameter(prs,'UseGpu',true,@(x) validateattributes(x,...
+    {'logical'},{'scalar'}));
+
+if ~iscell(Vi), Vi = {Vi}; end % Make sure Vi is a cell array
+parse(prs,Vi,BlkSz,Fun,varargin{:});
+VSz = size(Vi{1});
+NumIn = length(Vi);
+cellfun(@(x) validateattributes(x,...
+    {'numeric'},{'3d','real','size',VSz}),Vi);
+Opts = prs.Results;
+if Opts.UseGpu && gpuDeviceCount == 0
+    error('No GPU available');
+end
+
+if isscalar(BlkSz)
+    BlkSz = repmat(BlkSz,1,3);
+end
+if isscalar(Opts.BorderSize)
+    Opts.BorderSize = repmat(Opts.BorderSize,1,3);
+end
+
+if Opts.UseGpu
+    sendDat = @gpuArray;
+    recvDat = @gather;
+else
+    sendDat = @(x) x;
+    recvDat = @(x) x;
+end
+
+% Add additional 1 to end of sizes if needed, to ensure they are 3 vectors
+VSz = [VSz ones(1,3-length(VSz))];
+BlkSz = [BlkSz ones(1,3-length(BlkSz))];
+Opts.BorderSize = [Opts.BorderSize ones(1,3-length(Opts.BorderSize))];
+
+[B,E,BB,BE] = BlockIndices(VSz,BlkSz,Opts.BorderSize);
+
+nblks = length(B{1})*length(B{2})*length(B{3});
+blkdims = [length(B{1}),length(B{2}),length(B{3})];
+varargout = cell(Opts.NumOut,1);
+if nblks == 1
+    % In case we only need 1 block, just pass the entire volume to Fun.
+    % This is more efficient than using indexing since that takes a long 
+    % time for large inputs.
+    if Opts.Verbose
+        TextProgressBar('Processing blocks: ');
+        TextProgressBar(1, 1);
+    end
+    BlkIn = cellfun(sendDat,Vi,'UniformOutput',0);
+    [varargout{:}] = Fun(BlkIn{:});
+    varargout = cellfun(recvDat,varargout,'UniformOutput',0);
+else
+    if Opts.Verbose, TextProgressBar('Processing blocks: '); end
+    for oi = 1:Opts.NumOut
+        varargout{oi} = zeros(VSz,class(Vi{1}));
+    end
+
+    BlkResDevice = cell(1,Opts.NumOut);
+
+    if Opts.Verbose, TextProgressBar(1, nblks); end
+    for blki = 1:nblks
+        if Opts.Verbose, TextProgressBar(blki, nblks); end
+
+        % Prepare new block for sending
+        BlkInHost = GetHostBlock(blki);
+        % Send new block to the GPu
+        BlkInDevice = cellfun(sendDat,BlkInHost,'UniformOutput',0);
+        % Start processing for new block
+        [BlkResDevice{:}] = Fun(BlkInDevice{:});
+        % Transfer current result back from the GPU
+        BlkResHost = cellfun(recvDat,BlkResDevice,'UniformOutput',0);
+        % Place previous block result in output volume
+        PutHostBlock(blki,BlkResHost);
+        % Clear device results
+        for bi = 1:Opts.NumOut
+            BlkResDevice{bi} = [];
+        end
+    end
+end
+if Opts.Verbose, fprintf(' Done\n'); end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%                             LOCAL FUNCTIONS
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+    function BlkHost = GetHostBlock(blki)
+        [i,j,k] = ind2sub(blkdims, blki);
+        BlkHost = cell(1,NumIn);
+        for ii = 1:NumIn
+            BlkHost{ii} = Vi{ii}(...
+                BB{1}(i):BE{1}(i),...
+                BB{2}(j):BE{2}(j),...
+                BB{3}(k):BE{3}(k));
+        end
+    end
+
+    function PutHostBlock(blki,BlkResHost)
+        [i,j,k] = ind2sub(blkdims, blki);
+        for ii = 1:Opts.NumOut
+            varargout{ii}(...
+                B{1}(i):E{1}(i),...
+                B{2}(j):E{2}(j),...
+                B{3}(k):E{3}(k)) = BlkResHost{ii}(...
+                (B{1}(i)-BB{1}(i)+1):(end-(BE{1}(i)-E{1}(i))),...
+                (B{2}(j)-BB{2}(j)+1):(end-(BE{2}(j)-E{2}(j))),...
+                (B{3}(k)-BB{3}(k)+1):(end-(BE{3}(k)-E{3}(k))));
+        end
+    end
+end

ShapeMeasures.m

+function [Cl,Cs,Cp] = ShapeMeasures(E1,E2,E3)
+% SHAPEMEASURES Structure tensor shape measures
+% [Cl,Cs,Cp] = ShapeMeasures(E1,E2,E3)
+%
+% NOTE: Only the first nargout shape measures are actually computed. Thus
+% it is not necessary to provide the E3 input if only Cl is wanted.
+%
+% INPUT
+% E1, E2, E3 - Structure tensor eigenvalues from smallest to largest
+%
+% OUTPUT
+% Cl - Structure tensor line measure (needs E1 + E2)
+% Cs - Structure tensor isotropic measure (needs E1 + E3)
+% Cp - Structure tensor plane measure (needs E1 + E2 + E3)
+%
+% Patrick M. Jensen, 2018, ROCKWOOL
+
+% Cl = (E2 - E1)./(E2 + eps);
+% if nargout > 1
+%     Cs = E1./(E3 + eps);
+% end
+% if nargout > 2
+%     Cp = (E1.*(E3 - E2))./(E2.*E3 + eps);
+% end
+E1r = 1./(E1 + eps);
+E2r = 1./(E2 + eps);
+E3r = 1./(E3 + eps);
+
+Cl = (E1r - E2r)./E1r;
+Cp = (E2r - E3r)./E1r;
+Cs = E3r./E1r;
\ No newline at end of file

StructureTensor3.m

+function [T11,T22,T33,T12,T13,T23] = StructureTensor3(V,sigma,rho,varargin)
+% STRUCTURETENSOR3 Compute the structure tensors for a 3D volume
+% [T11,T22,T33,T12,T13,T23] = StructureTensor3(V,sigma,rho)
+% [T11,T22,T33,T12,T13,T23] = StructureTensor3(___,Name,Value,...)
+%
+% INPUT
+% V     - Input 3D volume of size N x M x L
+% sigma - Std. dev. of convolution kernel for derivatives
+% rho   - Std. dev. of convolution kernel for structure tensor
+% 
+% KEY-VALUE PARAMETERS:
+% Verbose   - Print status information (default: false)
+% UseGpu    - Use the GPU. Must be logical or text with 'auto'. If value is
+%             'auto' the GPU is used if any is available (default: 'auto')
+% BlockSize - Split the input into blocks of size BlockSize. This option is
+%             ignored if UseGpu is false. (default: size(V))
+%
+% OUTPUT
+% T11, T22, T33, T12 T13, T23 - N x M x L arrays with the elements of the
+%                               structure tensor
+%
+% See also EIGREALSYMM3
+%
+% Patrick M. Jensen, 2017, ROCKWOOL
+
+% Handle inputs
+narginchk(3,inf);
+prs = inputParser;
+% Required
+addRequired(prs,'V',@(x) validateattributes(x,...
+    {'double','single','gpuArray'},{'3d','nonsparse','real'}));
+addRequired(prs,'sigma',@(x) validateattributes(x,...
+    {'numeric'},{'scalar','real','positive'}));
+addRequired(prs,'rho',@(x) validateattributes(x,...
+    {'numeric'},{'scalar','real','positive'}));
+% Key-value parameters
+addParameter(prs,'Verbose',false,@(x) validateattributes(x,...
+    {'logical'},{'scalar'}));
+addParameter(prs,'UseGpu','auto',@(x) validateattributes(x,...
+    {'logical','string','char'},{}));
+addParameter(prs,'BlockSize',size(V),@(x) validateattributes(x,...
+    {'numeric'},{'integer','vector','positive'}));
+parse(prs,V,sigma,rho,varargin{:});
+Opts = prs.Results;
+if isscalar(Opts.BlockSize)
+    Opts.BlockSize = repmat(Opts.BlockSize,1,3);
+elseif length(Opts.BlockSize) ~= 3
+    error('BlockSize must be a scalar or 3-vector');
+end
+switch validatestring(string(Opts.UseGpu),{'false','true','auto'})
+    case 'false', Opts.UseGpu = false;
+    case 'true', Opts.USeGpu = true;
+    case 'auto', Opts.UseGpu = gpuDeviceCount > 0;
+end
+if Opts.UseGpu && gpuDeviceCount == 0
+    error('No GPU available');
+end
+
+[T11,T22,T33,T12,T13,T23] = Gpu3dBlockProc(V,Opts.BlockSize,...
+    @(V) CalcStructureTensor3(V,sigma,rho),'Verbose',Opts.Verbose,...
+    'UseGpu',Opts.UseGpu,'BorderSize',ceil(2*max(sigma,rho)) + 1,...
+    'NumOut',6);
+end

StructureTensorEig3.m

+function [E,X,Y,Z] = StructureTensorEig3(V,sigma,rho,varargin)
+% STRUCTURETENSOREIG3 Structure tensor eigenvalues and -vectors
+% [E,X,Y,Z] = StructureTensorEig3(V,sigma,rho)
+% [E,X,Y,Z] = StructureTensorEig3(___,Name,Value)
+% 
+% INPUT
+% V     - Input 3D volume of size N x M x L
+% sigma - Std. dev. of convolution kernel for derivatives
+% rho   - Std. dev. of convolution kernel for structure tensor
+%
+% KEY-VALUE PARAMETERS
+% NumEigVals - Return N smallest eigenvalues (default: 3)
+% NumEigVecs - Return eigenvectors for N smallest eigenvalues (default: 3)
+% Verbose    - Print status information (default: false)
+% UseGpu     - Use the GPU. Must be logical or text with 'auto'. If value
+%              is 'auto' the GPU is used if any is available 
+%              (default: 'auto')
+% BlockSize  - Block size for block processing (default: size(V))
+%
+% OUTPUT
+% E - NumEigVals x 1 cell with N x M x L arrays with eigenvalues
+% X - NumEigVecs x 1 cell with N x M x L array with the x-components
+% Y - NumEigVecs x 1 cell with N x M x L array with the y-components
+% Z - NumEigVecs x 1 cell with N x M x L array with the z-components
+%
+% See also STRUCTURETENSOR3 EIGREALSYMM3 GPU3DBLOCKPROC
+%
+% Patrick M. Jensen, 2018, ROCKWOOL
+
+narginchk(3,inf);
+prs = inputParser;
+% Required
+addRequired(prs,'V',@(x) validateattributes(x,...
+    {'double','single','gpuArray'},{'3d','nonsparse','real'}));
+addRequired(prs,'sigma',@(x) validateattributes(x,...
+    {'numeric'},{'scalar','real','positive'}));
+addRequired(prs,'rho',@(x) validateattributes(x,...
+    {'numeric'},{'scalar','real','positive'}));
+% Key-value parameters
+addParameter(prs,'Verbose',false,@(x) validateattributes(x,...
+    {'logical'},{'scalar'}));
+addParameter(prs,'UseGpu','auto',@(x) validateattributes(x,...
+    {'logical','string','char'},{}));
+addParameter(prs,'BlockSize',size(V),@(x) validateattributes(x,...
+    {'numeric'},{'integer','vector','positive'}));
+addParameter(prs,'NumEigVals',3,@(x) validateattributes(x,...
+    {'numeric'},{'integer','nonnegative','<=',3,'scalar'}));
+addParameter(prs,'NumEigVecs',3,@(x) validateattributes(x,...
+    {'numeric'},{'integer','nonnegative','<=',3,'scalar'}));
+parse(prs,V,sigma,rho,varargin{:});
+Opts = prs.Results;
+if isscalar(Opts.BlockSize)
+    Opts.BlockSize = repmat(Opts.BlockSize,1,3);
+elseif length(Opts.BlockSize) ~= 3
+    error('BlockSize must be a scalar or 3-vector');
+end
+switch validatestring(string(Opts.UseGpu),{'false','true','auto'})
+    case 'false', Opts.UseGpu = false;
+    case 'true', Opts.USeGpu = true;
+    case 'auto', Opts.UseGpu = gpuDeviceCount > 0;
+end
+if Opts.UseGpu && gpuDeviceCount == 0
+    error('No GPU available');
+end
+
+nval = Opts.NumEigVals;
+nvec = Opts.NumEigVecs;
+
+border = ceil(2*max(sigma,rho)) + 1;
+nout = nval+3*nvec;
+Res = cell(1,nout);
+[Res{:}] = Gpu3dBlockProc({V},Opts.BlockSize,...
+    @(V) ProcessBlock(V,sigma,rho,nval,nvec),...
+    'Verbose',Opts.Verbose,'NumOut',nout,'BorderSize',border,...
+    'UseGpu',Opts.UseGpu);
+E = cell(1,nval);
+X = cell(1,nvec);
+Y = cell(1,nvec);
+Z = cell(1,nvec);
+% Gather up results
+for i = 1:nval
+    E{i} = Res{i};
+end
+for i = 1:nvec
+    X{i} = Res{nval+(i-1)*3+1};
+    Y{i} = Res{nval+(i-1)*3+2};
+    Z{i} = Res{nval+(i-1)*3+3};
+end
+end
+
+function varargout = ProcessBlock(VBlock,sigma,rho,nval,nvec)
+varargout = cell(1,nval+3*nvec);
+[T11,T22,T33,T12,T13,T23] = CalcStructureTensor3(VBlock,sigma,rho);
+[varargout{:}] = CalcEigRealSymm3(T11,T22,T33,T12,T13,T23,nval,nvec);
+end
\ No newline at end of file

TextProgressBar.m

+function TextProgressBar(varargin)
+% TEXTPROGRESSBAR - Display a text progress bar
+% TextProgressBar(Init)
+% TextProgressBar(Crnt, Max)
+%
+% INPUT
+% Init - Initial part of line.
+% Crnt - Current value out of Max.
+% Max  - The maximum value.
+%
+% Patrick M. Jensen, 2018, ROCKWOOL
+
+% Modified version of: 
+% textprogressbar, 
+% https://se.mathworks.com/matlabcentral/fileexchange/28067-text-progress-bar
+
+persistent strCR;           %   Carriage return pesistent variable
+
+narginchk(1,2);
+
+% Vizualization parameters
+strPercentageLength = 10;   %   Length of percentage string (must be >5)
+strDotsMaximum      = 10;   %   The total number of dots in a progress bar
+
+% Main 
+Crnt = varargin{1};
+if ischar(Crnt)
+    % Progress bar - initialization
+    fprintf('%s',Crnt);
+    strCR = -1;
+elseif isnumeric(Crnt)
+    Max = varargin{2};
+    % Progress bar - normal progress
+    pct = floor(100 * Crnt / Max);
+    progressOut = [repmat(' ',1,floor(log10(Max)) - floor(log10(Crnt)))...
+        num2str(Crnt) '/' num2str(Max) ' | '];
+    percentageOut = [num2str(pct) '%%'];
+    percentageOut = [percentageOut...
+        repmat(' ',1,strPercentageLength-length(percentageOut)-1)];
+    nDots = floor(pct/100*strDotsMaximum);
+    dotOut = ['[' repmat('.',1,nDots) repmat(' ',1,strDotsMaximum-nDots) ']'];
+    strOut = [progressOut percentageOut dotOut];
+    
+    % Print it on the screen
+    if strCR == -1
+        % Don't do carriage return during first run
+        fprintf(strOut);
+    else
+        % Do it during all the other runs
+        fprintf([strCR strOut]);
+    end
+    
+    % Update carriage return
+    strCR = repmat('\b',1,length(strOut)-1);
+    
+else
+    % Any other unexpected input
+    error('Unsupported argument type');
+end
+
+end
\ No newline at end of file

VecOrientation.m

+function [Phi,Theta] = VecOrientation(X,Y,Z)
+% VECORIENTATION Vector azimuth (phi) and elevation (theta)
+% [Phi,Theta] = VecOrientation(X,Y,Z)
+%
+% Patrick M. Jensen, 2018, ROCKWOOL
+
+FlipMask = 1 - 2*(Z < 0);
+X = FlipMask.*X;
+Y = FlipMask.*Y;
+Z = FlipMask.*Z;
+
+Phi = atan2(-Y, X);
+Theta = atan2(Z, sqrt(X.^2 + Y.^2));
\ No newline at end of file

icosahedron_color.m

+function [R,G,B] = icosahedron_color(X,Y,Z)
+% Coloring of fiber orientation. 
+% 
+% Input parameters
+% v: 4D array representing the orientations in all voxels of the volume
+% mask: binary mask for volumetric data
+% saveDest: directory for saving the computed image data
+% slices: 1D array specifying the slices of which image data should be returned. 
+% 
+% S_display: colored image data for selected slices (as specified by the
+%            variable slices) 
+% 
+% Copyright 2017 MIT License
+% Vedrana Andersen Dahl (vand@dtu.dk) and Camilla H. Trinderup (ctri@dtu.dk)
+% Technical University of Denmark
+% 
+% Date created: 
+% 27/4-2017
+
+% dimension of volume
+[m,n,d,~] = size(X);
+
+% direction knot distances
+centers = cast(isosahedron_points(),'like',X);
+colors = hsv(6)';
+
+% compute coloring for all slices 
+S = zeros(m,n,d,3,'like',X);
+
+F = [X(:)'; Y(:)'; Z(:)'];
+distances = zeros(6, size(F,2),'like',F);
+for k = 1:6
+    distances(k,:) = undirected_distance(centers(k,:),F);
+end
+w = 1-distances/acos((1+sqrt(5))/(5+sqrt(5))); % weights, 1 for distance 0, 0 for maxdistance:
+w(w<0) = 0; % should not happen
+[w,i] = sort(w,'descend');
+sw = sum(w(1:3,:));
+sw(sw==0)=1; % for normalization so that weights sum to 1
+vc = repmat(w(1,:)./sw,[3,1]).*colors(:,i(1,:))...
+    + repmat(w(2,:)./sw,[3,1]).*colors(:,i(2,:))...
+    + repmat(w(3,:)./sw,[3,1]).*colors(:,i(3,:));
+
+% convert to rgb image
+R = reshape(vc(1,:),m,n,d);
+G = reshape(vc(2,:),m,n,d);
+B = reshape(vc(3,:),m,n,d);
+end
+
+function distances = undirected_distance(d1,data)
+arg = min(1, max(-1, d1*data)); % Truncate to keep acos real
+distances = real(min(acos(arg),acos(-arg)));
+
+end
+
+function S = isosahedron_points
+% Isosahedron points on a unit sphere
+
+phi=(1+sqrt(5))/2;
+
+S = [0 -1 -phi
+    0 -1 phi
+    -1 -phi 0
+    -1 phi 0
+    -phi 0 -1
+    phi 0 -1];
+
+S = S/sqrt(1+phi^2);
+
+end
\ No newline at end of file