Code explanation

Interested in generating hybrid FODs with own diffusion MRI and microscopy data? Here we provide details of the method alongside important snippets of code. The full code is made openly at https://git.fmrib.ox.ac.uk/srq306/hybrid_microscopy_dmri.

Diffusion MRI orientation

Diffusion MRI (dMRI) orientations were estimated using the Ball and stick model in FSL. The signal decay of the diffusion was modelled by an isotropic compartment (Ball) and multiple restricted compartments for fibre orientation (Stick).

dMRI was analysed to estimate ≤ 3 fibre populations with 50 orientation estimates (samples) per population.

_images/fig16.jpg

For each fibre population <i>, BAS model outputs three files:
merged_th<i>samples:
merged_ph<i>samples:
merged_f<i> samples:
where th and ph describe the 3D fibre orientation and f is the signal fraction (sometimes called voume fraction, though this ignores T2 effects)

Fibre populations with signal fractions less than 0.05 were removed. Details of ball and stick model output can be found here

Microscopy orientation

First, the microscopy orientation was registered to the dMRI space with TIRL and FSL FLIRT. The voxel coordinate for each microscopic pixel was then described by x, y, z coordinates in the diffusion mri space (c1,c2,c3). The fiber orientation was represented by a 3D vector (v1,v2,v3).

_images/fig1.jpg

Second, to facilitate a comparison between the 3D dMRI fibre orientation and 2D microscopy orientation, the dMRI orientation (a) was projected onto the 2D microscopy plane (a1) and onto the normal vector of the plane (a2). Next, the angle difference between the projected BAS orientation and the microscopy orietnation was calculated.

_images/fig2.jpg

function out = project(v,n)
%Ensure norm is unit
n = n./vecnorm(n,2,2);
% Find vector projection
a1 = sum(v.*n,2).*n;
% i.e. the vector projected onto the plane
a2 = v-a1;
out.a1 = a1;
out.a2 = a2;
end

Third, the BAS sample with the smallest angle to the microscopy orientation on the microscopy plane was selected. The through plane angle of the diffusion sample was used for the hybrid orientation reconstruction.

_images/fig3.jpg

 % Record smallest angle between vector and any dyad sample
[~,indd] = max(cosangsqrd,[],2,'omitnan');
linearind = sub2ind(size(cosangsqrd),1:size(cosangsqrd,1),indd'); %'
angg = acos(sqrt(cosangsqrd(linearind)));


a1 = reshape(a1,[],3);
a2 = reshape(a2,[],3);

% Output dyad sample most closely associated with each micro orientation
selected.a1 = a1(linearind,:);
selected.a2 = a2(linearind,:);
selected.ang = angg;
selected.ind = indd;

Hybrid orientation

We merged the in-plane orientation from the microscopy data with the through-plane orientation of the best matched diffusion sample to create a 3D hybrid orientation.

_images/fig4.jpg

micro.vect3D = tmp.inplane.*vecnorm(tmp.a2,2,2)+tmp.a1;
micro.vect3D = reshape(micro.vect3D,size(micro.inplane));

We then combine the hybrid fibre orientations over a local neighbourhood with the size of the desired voxel to create a hybrid fibre orientation distribution (FOD). For each voxel, the orientation was compared to a 3D vector set (256 directions evenly distributed across a sphere) to populate a 3D frequency histogram.

% For each voxel in hr space, extract fibre orientations, compare
% to directions in 'vectors' and populate frequency histogram
VV = unique(voxind);
VV(isnan(VV)) = [];
for w = 1:numel(VV)
    v = VV(w);
    if roimask_us(v)==1
        ind = voxind==v;
        if sum(ind)>0
            out = hist_sphere(micro.vect3D(ind,:),vectors);
            count(v,:) = count(v,:)+out.count;
        end
    end
    if mod(w,500)==0, disp([num2str(w) '/' num2str(numel(VV))]); end
 end

Spherical harmonics were then fitted to to the frequency histogram to generate our hybrid FODs.

% Fit SH coeffs to frequency matrix
if mrtrixflag
    vectors(1,:) = -vectors(1,:);
    disp('LR flip for mrtrix')
end

SHmat = SH_transform(vectors,8); % get the spherical harmonics basis
SHmat_pinv = pinv(SHmat);
Ncoeffs = size(SHmat,2);

% Normalise hsitogram by number of pixels in each voxel
countn = count./sum(count,2);
coeffs = SHmat_pinv*countn';
SH_3D = reshape(coeffs',s1,s2,s3,Ncoeffs);
count_3D = reshape(count,s1,s2,s3,256);