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.

.. sidebar:: Note

    The code provided in this page is not written as a one-line command for general usage. Instead, it serves as a walk-through guide on how to generate hybrid FODs. This guide is particularly helpful for adapting the code to suit your own study. If you're interested in understanding the method behind it, feel free to ignore the code block.



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. 

.. image:: ../html/_images/fig16.jpg
  :align: center
  :width: 300px

|

| 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 <https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/FDT/UserGuide>`_

.. sidebar:: Note

    Other diffusion models providing fibre orientation estimates can also be used.



Microscopy orientation
--------
First, the microscopy orientation was registered to the dMRI space with `TIRL <https://git.fmrib.ox.ac.uk/ihuszar/tirl>`_ and FSL `FLIRT <https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/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)**.

.. image:: ../html/_images/fig1.jpg
  :align: center
  :width: 400px
|

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.

.. image:: ../html/_images/fig2.jpg
  :align: center
  :width: 300px
|

.. code-block:: matlab


        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.

.. image:: ../html/_images/fig3.jpg
  :align: center
  :width: 300px
|

.. code-block:: matlab
 
         % 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.

.. image:: ../html/_images/fig4.jpg
  :align: center
  :width: 300px
|

.. code-block:: matlab
        
       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.

.. code-block:: matlab

        % 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.

.. code-block:: matlab

        % 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);