feComputeCanonicalDiffusion

 Calculate a tensor for forward modeling at each point of each fiber.

  Q = feComputeCanonicalDiffusion(fibers,dParms)

 INPUTS:
   fibers  - is a cell array of x,y,z coordinates for each node in a fiber.
   dParams - is a a vector of axial and radial diffusivity to use as
              parameters of the canonical tensors used to compute the 
              diffusion at the node.
 OUTPUT:
   Q       - is a cell array of the same length as the number of fibers.
             Each Q{ii} contains a matrix of (numNodes x 9) tensors.

 To put the tensor into the quadratic form, use T = reshape(T,3,3);
 eigs(T) calculates the axial diffusivity (largest) and so forth.

 EXAMPLE:
  d_ad = 1.5; d_rd = 0.3;
  dParms(1) = d_ad; dParms(2) = d_rd; dParms(3) = d_rd;
  fgImg.Q = feComputeCanonicalDiffusion(fibers,dParms)

 See also: feConnectomeInit.m v_lifeExample.m

 Copyright (2013-2014), Franco Pestilli, Stanford University, pestillifranco@gmail.com.

0001 function Q = feComputeCanonicalDiffusion(fibers,dParms)
0002 % Calculate a tensor for forward modeling at each point of each fiber.
0003 %
0004 %  Q = feComputeCanonicalDiffusion(fibers,dParms)
0005 %
0006 % INPUTS:
0007 %   fibers  - is a cell array of x,y,z coordinates for each node in a fiber.
0008 %   dParams - is a a vector of axial and radial diffusivity to use as
0009 %              parameters of the canonical tensors used to compute the
0010 %              diffusion at the node.
0011 % OUTPUT:
0012 %   Q       - is a cell array of the same length as the number of fibers.
0013 %             Each Q{ii} contains a matrix of (numNodes x 9) tensors.
0014 %
0015 % To put the tensor into the quadratic form, use T = reshape(T,3,3);
0016 % eigs(T) calculates the axial diffusivity (largest) and so forth.
0017 %
0018 % EXAMPLE:
0019 %  d_ad = 1.5; d_rd = 0.3;
0020 %  dParms(1) = d_ad; dParms(2) = d_rd; dParms(3) = d_rd;
0021 %  fgImg.Q = feComputeCanonicalDiffusion(fibers,dParms)
0022 %
0023 % See also: feConnectomeInit.m v_lifeExample.m
0024 %
0025 % Copyright (2013-2014), Franco Pestilli, Stanford University, pestillifranco@gmail.com.
0026 
0027 % Preallocate
0028 nFibers = length(fibers); % The number of Fibers.
0029 Q       = cell(1,nFibers); % Memory for the tensors of each fiber.
0030 D       = diag(dParms);    % The diagonal form of the Tensors' model parameters.
0031 
0032 feOpenLocalCluster
0033 parfor ii = 1:nFibers
0034  % Compute the diffusion gradient at each node of the fiber.
0035  fiberGradient = gradient(fibers{ii});
0036  
0037  % Number of nodes fro this fiber
0038  numNodes = size(fibers{ii},2);
0039  
0040  % preallocated memory for the vector representation of tensors.
0041  T = zeros(numNodes,9);
0042  
0043  for jj = 1:numNodes
0044   % Rotate the tensor toward the gradient of the fiber.
0045   %
0046   % Calculate a rotation matrix for the tensor so that points in the fiberGradient
0047   % direction and has two perpendicular directions (U)
0048   % Leaving the 3 outputs for this function is the fastest use of it.
0049   [Rot,~, ~] = svd(fiberGradient(:,jj)); % Compute the eigen vectors of the gradient.
0050   
0051   % Create the quadratic form of the tensor.
0052   %
0053   % The principal eigenvector is in the same direction of the
0054   % fiberGradient. The direction of the other two are scaled by dParms.
0055   % Human friendly version fo the code:
0056   % tensor = Rot*D*Rot'; % tensor for the current node, 3x3 matrix.
0057   % T(jj,:) = reshape(tensor,1,9); % reshaped as a vector 1,9 vector
0058   T(jj,:) = reshape(Rot*D*Rot',1,9);
0059  end
0060  % T is a matrix; each row is a 1,9 version of the tensor.
0061  Q{ii} = T;
0062 end
0063 
0064 return
0065