LIE ANALYSIS TUTORIAL

Enhancing Diffusion's on 3D Images

Note this notebook is a work in progress, and does not make optimal use of the LieAnalysis package.

LeftInvariantDerivatives3D[Obj3DPositionOrientationData, derivativeIndex]	computes the A_{derivativeIndex} derivatives from Obj3DPositionOrientationData using a finite differences method
LeftInvariantDerivatives3D[Obj3DPositionOrientationData, {σs,σo}, derivativesIndex]	using Gaussian derivatives specified by the scales σs and σo
OrientationScoreTensor3D[Obj3DPositionOrientationData, {derivativeIndex..}]	construct tensor {A_{derivativeIndex}..} from Obj3DPositionOrientationData using finite differences method
OrientationScoreTensor3D[Obj3DPositionOrientationData, {σs,σo}, {derivativeIndex..}]	using Gaussian derivatives specified by the scales σs and σo
OrientationScoreGaugeFrames3D[Obj3DPositionOrientationData,{σs,σo}]	computes a data adaptive frame from Obj3DPositionOrientationData
OrientationScoreTransform3D[data]	computes an orientation score (Obj3DPositionOrientationData) from data

Functions that can be used for enhancing diffusion's on 3-dimensional images.

Load the package

In[1]:=

Helper Functions

You should run this section, but you might not want to read this section.

Functions to compute structureness

A lot of the functions bellow are included in the LieAnalysis package internally, and should be removed in future releases of this Tutorial.

RotateZY[n_] := RotateZY@@Transpose[n,Join[Range[2, ArrayDepth[n]], {1}]];
RotateZY[x_,y_,z_] := RotateZYC[x,y,z];
RotateZYC = Compile[{{x,_Real},{y,_Real},{z,_Real}},Module[{β,γ,ϵ=10^-8},
    If[z<1-ϵ,
        If[z>-1+ϵ,
            {β,γ}={ArcCos[z],ArcTan[x,y]};,
            {β,γ}={π,0};
        ];,
        {β,γ}={0,0};
    ];
    {
        {Cos[β] Cos[γ],-Sin[γ],Cos[γ] Sin[β]},
        {Cos[β] Sin[γ],Cos[γ],Sin[β] Sin[γ]},
        {-Sin[β],0,Cos[β]}
    }],
    RuntimeAttributesListable
];

CartesianToSpherical[qvecs_] := CartesianToSphericalC[qvecs];
CartesianToSphericalC=Compile[{{vec,_Real,1}},
    Module[{qx,qy,qz,R,θ,ϕ},
        {qx,qy,qz} = vec; R = Sqrt[qx^2+qy^2+qz^2];
        {qx,qy,qz} = {qx,qy,qz}/R;
        If[qz<1-10^-8,
            If[qz>-1+10^-8,
                {θ,ϕ} = {ArcCos[qz/Sqrt[qx^2+qy^2+qz^2]],ArcTan[qx,qy]};,
                {θ,ϕ} = {π,0};
            ];,
            {θ,ϕ} = {0,0};
        ];
        {R,θ,ϕ}
    ],
    RuntimeAttributes->{Listable}
];

SHIndexj2l[j_] :=  Floor[Sqrt[j-1]];
SHIndexlm2j[l_,m_] := Round[(l^2+l+1)+m];
SHIndexj2m[j_] := Round[j-1-(SHIndexj2l[j])^2-SHIndexj2l[j]];

DSHTCreateM[orientations_,lMax_Integer] := Module[{Mmat,R,θ,ϕ},
    {R,θ,ϕ} = CartesianToSpherical[orientations];
    Mmat = Table[SphericalHarmonicY[SHIndexj2l[j],SHIndexj2m[j],θ,ϕ ],{j,1,(lMax+1)^2}];
    Mmat
]

PseudoInverseiDSHT[orientations_,Lmax_Integer,values_] := values.Transpose[PseudoInverseiDSHT[orientations,Lmax]];
PseudoInverseiDSHT[orientations_,Lmax_Integer] :=  Module[{Mmat=DSHTCreateM[orientations,Lmax]},
    PseudoInverse[Transpose[Mmat]]
]

SHAngularDiffusion[s_List,sigma_] := s.Transpose[SHAngularDiffusion[SHIndexj2l[Dimensions[s][[-1]]],sigma]];
SHAngularDiffusion[lMax_Integer,sigma_] := Module[{kt = 1/2 sigma^2},
    (*If[E^(-lMax(lMax+1)kt)>0.01,Message[SHAngularDiffusion::nesh]];*)
    DiagonalMatrix[Table[E^(-SHIndexj2l[jp](SHIndexj2l[jp]+1)kt),{jp,1,SHIndexlm2j[lMax,lMax]}]]
]

(*Create GaussianDerivative using gaussian filter*)
MGaussianDerivative[data_,σGaus_,der_] := GaussianFilter[data,{5σGaus,σGaus},der,Method->"Gaussian","Standardization"->False];

(*Apply a gaussian filter with size σGaus at each orientation. der is the derivative you want to take with {0,0,0} for blurring.*)
GaussianFilterAtOrientations[OS_,σGaus_,der_] := Transpose[Table[MGaussianDerivative[OS[[All,All,All,i]],σGaus,der],{i,1,Length[OS[[1,1,1]]]}],{4,1,2,3}];

(*Calculate the spatial gradient at each orientation by applygin GaussianFilterAtOrientations 3 times*)
SpatialGradientAtOrientations[OS_,σGaus_] := Map[GaussianFilterAtOrientations[OS,σGaus,#]&,{{1,0,0},{0,1,0},{0,0,1}}];

(*Calculate the elements of the spatial Hessian by applying gaussianfilter 6 times*)
SpatialHessianAtOrientations[OS_,σGaus_] := Module[{dxx,dxy,dxz,dyy,dyz,dzz},
    dxx = GaussianFilterAtOrientations[OS,σGaus,{2,0,0}];
    dxy = GaussianFilterAtOrientations[OS,σGaus,{1,1,0}];
    dxz = GaussianFilterAtOrientations[OS,σGaus,{1,0,1}];

    dyy = GaussianFilterAtOrientations[OS,σGaus,{0,2,0}];
    dyz = GaussianFilterAtOrientations[OS,σGaus,{0,1,1}];

    dzz = GaussianFilterAtOrientations[OS,σGaus,{0,0,2}];
    {dxx,dxy,dxz,dyy,dyz,dzz}
]

HRegularizeAndPrepareSpatialDerivativesXYZ[s_, lMax_, sigmaSpatial_, sigmaAngle_]  := Module[{sAngularlyBlurred, sV, pack = Developer`ToPackedArray, sGradientXYZ, sHessianXYZ}, 
    (*Perform angular blurring*)
    sAngularlyBlurred = SHAngularDiffusion[s, sigmaAngle];

    (*Perform spatial blurring*)
    sV = GaussianFilterAtOrientations[sAngularlyBlurred, sigmaSpatial, {0, 0, 0}];

    (*Determine spatial derivatives in global coordinates x, y, z*)
    sGradientXYZ  = pack[SpatialGradientAtOrientations[sAngularlyBlurred, sigmaSpatial]];
    sHessianXYZ = pack[SpatialHessianAtOrientations[sAngularlyBlurred, sigmaSpatial]];

    {sV, sGradientXYZ, sHessianXYZ}
]

Options[CompiledSHTable] = {"PhiSymmetric" -> False};
CompiledSHTable[lMax_Integer, OptionsPattern[]] :=  Compile[{th, ph}, 
    Evaluate[
        If[OptionValue["PhiSymmetric"], 
            N[ExpToTrig[Table[If[SHIndexj2m[j] == 0, SphericalHarmonicY[SHIndexj2l[j], SHIndexj2m[j], th, ph], 0], {j, 1, (lMax+1)^2}]]], 
            N[ExpToTrig[Table[SphericalHarmonicY[SHIndexj2l[j], SHIndexj2m[j], th, ph], {j, 1, (lMax+1)^2}]]]
        ]
    ], 
    RuntimeAttributes -> Listable
];

CompiledSHA1A2A3[lMax_Integer, selection_:All] := 
    With[{csht = CompiledSHTable[lMax], sel = selection/.{"A1"->1, "A2"->2, "A3"->3, All->{1, 2, 3}}}, 
        (*Check option input*)
        If[(!ContainsOnly[sel, {1, 2, 3}])||Depth[sel] != 2, Abort[]];
    
        (*Compile function*)
        Compile[{{Rn1, _Real, 2}, {sGradR3, _Complex, 2}}, 
            Module[{n1 = Rn1.{0, 0, 1}, th, ph, A123}, 
                (*Calculate spherical angles of n1*)
                {th, ph} = CartesianToSphericalC[n1][[2;;3]];
        
                (*Calculate first order spatial derivatives*)
                A123 = Transpose[Rn1[[All, sel]]].(sGradR3.csht[th, ph])
            ], 
        RuntimeAttributes->Listable, CompilationOptions->{"InlineCompiledFunctions"->True, "InlineExternalDefinitions"->True}]
    ];

HHessianXYZToHessianA1A2A3CMatrix = Compile[{{sel, _Integer, 1}, {R, _Real, 2}}, 
    Module[{R11 = R[[1, 1]], R12 = R[[1, 2]], R13 = R[[1, 3]], R21 = R[[2, 1]], R22 = R[[2, 2]], R23 = R[[2, 3]], R31 = R[[3, 1]], R32 = R[[3, 2]], R33 = R[[3, 3]]}, 
        {
            {R11^2, 2 R11 R21, 2 R11 R31, R21^2, 2 R21 R31, R31^2}, 
            {R11 R12, R12 R21+R11 R22, R12 R31+R11 R32, R21 R22, R22 R31+R21 R32, R31 R32}, 
            {R11 R13, R13 R21+R11 R23, R13 R31+R11 R33, R21 R23, R23 R31+R21 R33, R31 R33}, 
            {R12^2, 2 R12 R22, 2 R12 R32, R22^2, 2 R22 R32, R32^2}, 
            {R12 R13, R13 R22+R12 R23, R13 R32+R12 R33, R22 R23, R23 R32+R22 R33, R32 R33}, 
            {R13^2, 2 R13 R23, 2 R13 R33, R23^2, 2 R23 R33, R33^2}
        }[[sel]]
    ], 
    RuntimeAttributes->{Listable}
];

CompiledSHSpatialHessian[lMax_Integer, selection_:All] := With[{csht = CompiledSHTable[lMax], sel = selection/.{"A1A1"->1, "A1A2"->2, "A1A3"->3, "A2A2"->4, "A2A3"->5, "A3A3"->6, All->{1, 2, 3, 4, 5, 6}}},     
    (*Check option input*)
    If[(!ContainsOnly[sel, {1, 2, 3, 4, 5, 6}])||Depth[sel] != 2, Abort[]];

    (*Compile function*)
    Compile[{{Rn1, _Real, 2}, {sHessianxyz, _Complex, 2}}, 
        Module[{n1 = Rn1.{0, 0, 1}, th, ph, A11A12A13A22A23A33}, 
            (*Calculate spherical angles of n1*)
            {th, ph} = CartesianToSphericalC[n1][[2;;3]];
    
            (*Calculate second order spatial derivatives*)
            A11A12A13A22A23A33 = HHessianXYZToHessianA1A2A3CMatrix[sel, Rn1].(sHessianxyz.csht[th, ph])
        ], 
        RuntimeAttributes->Listable, CompilationOptions->{"InlineCompiledFunctions"->True, "InlineExternalDefinitions"->True}
    ]
];

Ytilde[l_, m_, p_,θ_,ϕ_] := ^(- p ϕ) SphericalHarmonicY[l, m+p,θ,ϕ];

k[l_, m_] := -2 (l^2 + l - m^2 );
j[l_, m_] := Sqrt[(l + m + 2) (l - m - 1) (l + m + 1) (l - m)];
h[l_, m_] := 2 m Sqrt[(l - m) (l + m)];
g[l_, m_] := Sqrt[(l - m - 2) (l - m - 1) (l - m) (l + m + 1)];
a[l_, m_] := Sqrt[(l + m) (l + m - 1)];
b[l_, m_] := Sqrt[(l - m) (l + m + 1)];

shA4[l_,m_,θ_,ϕ_] :=  * 1/2 Sqrt[(2 l + 1)/(2 l - 1)] (Ytilde[l - 1, m, - 1 ,θ,ϕ] a[l, m] + Ytilde[l - 1, m, + 1,θ,ϕ] a[l, -m]);

shA5[l_,m_,θ_,ϕ_]:= 1/2 (b[l, m] (Ytilde[l, m, + 1,θ,ϕ]) -   b[l, -m] (Ytilde[l, m, - 1,θ,ϕ]) );

shA44[l_,m_,θ_,ϕ_]:= -l (l + 1) SphericalHarmonicY[l, m,θ,ϕ] - shA55[l, m ,θ,ϕ];

shA45[l_,m_,θ_,ϕ_]:= (Sqrt[2 l + 1])/ (4 Sqrt[2 l - 1]) (h[l, m] SphericalHarmonicY[l - 1, m,θ,ϕ] - g[l, -m] Ytilde[l - 1, m, -2,θ,ϕ] + g[l, m] Ytilde[l - 1, m, +2,θ,ϕ]);

shA55[l_,m_,θ_,ϕ_] := 1/4  (k[l, m]  Ytilde[l, m,0,θ,ϕ] + j[l, -m] Ytilde[l, m, -2,θ,ϕ] +  j[l, m] Ytilde[l, m, +2,θ,ϕ] );

SHA4[lMax_Integer, θ_, ϕ_] := Module[{M,beta,gamma},
    Table[shA4[SHIndexj2l[jp],SHIndexj2m[jp], θ, ϕ],{jp,1,SHIndexlm2j[lMax,lMax]}]
]

SHA5[lMax_Integer, θ_, ϕ_] := Module[{M,beta,gamma},
    Table[shA5[SHIndexj2l[jp],SHIndexj2m[jp], θ, ϕ],{jp,1,SHIndexlm2j[lMax,lMax]}]
]
SHA44[lMax_Integer, θ_, ϕ_] := Module[{M,beta,gamma},
    Table[shA44[SHIndexj2l[jp],SHIndexj2m[jp], θ, ϕ],{jp,1,SHIndexlm2j[lMax,lMax]}]
]

SHA55[lMax_Integer, θ_, ϕ_] := Module[{M,beta,gamma},
    Table[shA55[SHIndexj2l[jp],SHIndexj2m[jp], θ, ϕ],{jp,1,SHIndexlm2j[lMax,lMax]}]
]

SHA45[lMax_Integer, θ_, ϕ_] := Module[{M,beta,gamma},
    Table[shA45[SHIndexj2l[jp],SHIndexj2m[jp], θ, ϕ],{jp,1,SHIndexlm2j[lMax,lMax]}]
]

CompiledSHA4[lMax_Integer] := Compile[{th,ph},
    Evaluate[SHA4[lMax, th, ph]],
    RuntimeAttributesListable
];

CompiledSHA5[lMax_Integer] := Compile[{th,ph},
    Evaluate[SHA5[lMax, th, ph]],
    RuntimeAttributesListable
];

CompiledSHA44[lMax_Integer] := Compile[{th,ph},
    Evaluate[SHA44[lMax, th, ph]],
    RuntimeAttributesListable
];

CompiledSHA45[lMax_Integer] := Compile[{th,ph},
    Evaluate[SHA45[lMax, th, ph]],
    RuntimeAttributesListable
];

CompiledSHA55[lMax_Integer] := Compile[{th,ph},
    Evaluate[SHA55[lMax, th, ph]],
    RuntimeAttributesListable
];

SHRotateZYZ[gamma_,beta_,alpha_,s_List] := s.Transpose[SHRotateZYZ[gamma,beta,alpha,SHIndexj2l[Dimensions[s][[-1]]]]];
SHRotateZYZ[gamma_,beta_,alpha_,lMax_Integer] := Module[{M},
    M = ConstantArray[0,{SHIndexlm2j[lMax,lMax],SHIndexlm2j[lMax,lMax]}];
    Table[Table[(M[[jp,SHIndexlm2j[SHIndexj2l[jp],mp]]]=WignerD[{SHIndexj2l[jp],SHIndexj2m[jp],mp},-gamma,-beta,-alpha]),{mp,-SHIndexj2l[jp],SHIndexj2l[jp]}],{jp,1,SHIndexlm2j[lMax,lMax]}];
    M
]

RToAlphaBetaGamma=Compile[{{R,_Real,2}},
    Module[{γ,β,α,cb},
        cb=Chop[R[[3,3]]];
        If[cb<1-10^-8,
            If[cb>-1+10^-8,
                β = ArcCos[cb];
                γ = ArcTan[R[[1,3]],R[[2,3]]];
                α = ArcTan[-R[[3,1]],R[[3,2]]];,
                γ = 0; β = π; α = ArcTan[R[[2,2]],R[[2,1]]];
            ],
            γ=β =0;α=ArcTan[R[[1,1]],R[[2,1]]];
        ];
        {γ,β,α}
    ],RuntimeAttributes->{Listable}
]

Out[49]=

SHRotateR[R_,s_List] := s.Transpose[SHRotateR[R,SHIndexj2l[Dimensions[s][[-1]]]]];
SHRotateR[R_,lMax_Integer] := Module[{M,j,l,m,beta,gamma,alpha},
    {gamma,beta,alpha}=RToAlphaBetaGamma[R];
SHRotateZYZ[gamma,beta,alpha,lMax]]

Unprotect[Rez,Rey,Rex,ex,ey,ez]
ez={0,0,1};
ey={0,1,0};
ex={1,0,0};
Rez[α_] := RotationMatrix[α, ez];
Rey[γ_] := RotationMatrix[γ, ey];
Rex[β_] := RotationMatrix[β, ex];

Protect[Rez,Rey,Rex,ex,ey,ez]

Out[52]=

Out[59]=

RotateZYAngles=Compile[{{β,_Real},{γ,_Real}},
    {
        {Cos[β] Cos[γ],-Sin[γ],Cos[γ] Sin[β]},
        {Cos[β] Sin[γ],Cos[γ],Sin[β] Sin[γ]},
        {-Sin[β],0,Cos[β]}
    },
    RuntimeAttributesListable
];

CompiledSHA4A5[lMax_Integer, selection_:All] := 
    With[{csh4 = CompiledSHA4[lMax], csh5 = CompiledSHA5[lMax], Sey = SHRotateR[Rey[-π/2], lMax], Ry = Rey[π/2], ez = {0, 0, 1}, sel = selection/.{"A4"->1, "A5"->2, All->{1, 2}}}, 
        (*Check option input*)
        If[(!ContainsOnly[sel, {1, 2}])||Depth[sel] != 2, Abort[]];
            
        (*Compile function*)
        Compile[{{Rn1, _Real, 2}, {s, _Complex, 1}}, 
            Module[{n1 = Rn1.ez, closeToEz, th, ph, A4, A5, RalphaTilde, Rnα0, n2}, 
                (*Calculate spherical angles*)
                {th, ph} = CartesianToSphericalC[n1][[2;;3]];
                Rnα0 = RotateZYAngles[th, ph];

                (*Calculate alpha difference RalphaTilde*)
                RalphaTilde = Transpose[Rnα0].Rn1;
        
                (*Calculate first order angular derivatives*)
                A4 = s.csh4[th, ph];
                A5 = s.csh5[th, ph];
                
                (*Rotate using RalphaTilde*)
                (RalphaTilde[[1;;2, 1;;2]].{A4, A5})[[sel]]
            ], 
        RuntimeAttributes->Listable, CompilationOptions->{"InlineCompiledFunctions"->True, "InlineExternalDefinitions"->True}]
    ];

CompiledSHA44A45A55[lMax_Integer, selection_:All] := 
    With[{csh44 = CompiledSHA44[lMax], csh45 = CompiledSHA45[lMax], csh55 = CompiledSHA55[lMax], Sey = SHRotateR[Rey[-π/2], lMax], Ry = Rey[π/2], ez = {0, 0, 1}, sel = selection/.{"A4A4"->1, "A4A5"->2, "A5A5"->4, All->{1, 2, 4}}}, 
        (*Check option input*)
        If[(!ContainsOnly[sel, {1, 2, 4}])||Depth[sel] != 2, Abort[]];
            
        (*Compile function*)
        Compile[{{Rn1, _Real, 2}, {s, _Complex, 1}}, 
            Module[{n1 = Rn1.ez, th, ph, A44, A45, A55, RalphaTilde, Rnα0}, 
                (*Calculate spherical angles*)
                {th, ph} = CartesianToSphericalC[n1][[2;;3]];
                Rnα0 = RotateZYAngles[th, ph];

                (*Calculate alpha difference RalphaTilde*)
                RalphaTilde = Transpose[Rnα0].Rn1;

                (*Calculate second order angular derivatives*)
                A44 = s.csh44[th, ph];
                A45 = s.csh45[th, ph];
                A55 = s.csh55[th, ph];

                (*Rotate using RalphaTilde*)
                Flatten[(RalphaTilde[[1;;2, 1;;2]]).{{A44, A45}, {A45, A55}}.(Transpose[RalphaTilde][[1;;2, 1;;2]]), 1][[sel]]
            ], 
        RuntimeAttributes->Listable, CompilationOptions->{"InlineCompiledFunctions"->True, "InlineExternalDefinitions"->True}]
    ];

CompiledSHMixedDerivativesFromGradientR3[lMax_Integer, selection_:All] := 
    With[{csh4 = CompiledSHA4[lMax], csh5 = CompiledSHA5[lMax], Sey = SHRotateR[Rey[-π/2], lMax], Ry = Rey[π/2], ez = {0, 0, 1}, sel = selection/.{"A1A4"->1, "A2A4"->2, "A3A4"->3, "A1A5"->4, "A2A5"->5, "A3A5"->6, All->{1, 2, 3, 4, 5, 6}}}, 
        (*Check option input*)
        If[(!ContainsOnly[sel, {1, 2, 3, 4, 5, 6}])||Depth[sel] != 2, Abort[]];
        
        (*Compile function*)
        Compile[{{Rn1, _Real, 2}, {sGradR3, _Complex, 2}}, 
            Module[{n1 = Rn1.ez, th, ph, A4, A5, RalphaTilde, Rnα0}, 
                (*Calculate spherical angles*)
                {th, ph} = CartesianToSphericalC[n1][[2;;3]];
                Rnα0 = RotateZYAngles[th, ph];

                (*Calculate alpha difference RalphaTilde*)
                RalphaTilde = Transpose[Rnα0].Rn1;
        
                (*Calculate first order angular derivatives*)
                A4 = sGradR3.csh4[th, ph];
                A5 = sGradR3.csh5[th, ph];
                
                (*Rotate using RalphaTilde*)
                Flatten[RalphaTilde[[1;;2, 1;;2]].{A4, A5}.(Rn1)][[sel]]
            ], 
        RuntimeAttributes->Listable, CompilationOptions->{"InlineCompiledFunctions"->True, "InlineExternalDefinitions"->True}]
    ];

HCalculateDerivatives::usage = "HCalculateDerivatives[s,\!\(\*SuperscriptBox[\(s\), \(∇\)]\),\!\(\*SuperscriptBox[\(s\), \(H\)]\), R, L , derivativeList] calculates the derivatives specified by derivativeList from s(x,·): the spherical harmonic coefficients of V(x,·), \!\(\*SuperscriptBox[\(s\), \(∇\)]\)(·,x,·): the spherical harmonic coefficients of \!\(\*SubscriptBox[\(∇\), SuperscriptBox[\(\), \(3\)]]\)V(x,·), \!\(\*SuperscriptBox[\(s\), \(H\)]\)(·,x,·): the spherical harmonic coefficients of \!\(\*SubscriptBox[\(H\), SuperscriptBox[\(R\), \(3\)]]\)V(x,·) and list of 3D rotations R\nInput should have the following dimensions [s]={\!\(\*SubscriptBox[\(N\), \(pos\)]\),j(L,L)}, [\!\(\*SuperscriptBox[\(s\), \(∇\)]\)]={3, \!\(\*SubscriptBox[\(N\), \(pos\)]\),j(L,L)}, [\!\(\*SuperscriptBox[\(s\), \(H\)]\)]={6,\!\(\*SubscriptBox[\(N\), \(pos\)]\),j(L,L)}, [R]={\!\(\*SubscriptBox[\(N\), \(pos\)]\),\!\(\*SubscriptBox[\(N\), \(o\)]\),3,3}, and derivativeList should be a list containing all or a subset of {'A1','A2','A3','A1A1','A1A2','A1A3','A2A2','A2A3','A3A3','A4','A5','A4A4','A4A5','A5A5','A1A4','A2A4','A3A4','A1A5','A2A5','A3A5'}";
HCalculateDerivatives[sV_, sGradientXYZ_, sHessianXYZ_, Rn_, lMax_, derList_]  := Module[{mixedDer, firstOrderSpatialDerivatives, secondOrderSpatialDerivatives, firstOrderAngularDerivatives, secondOrderAngularDerivatives, mixedDerivatives, A1A2A3, A4A5, A44A45A55, A11A12A13A22A23A33}, 
    firstOrderSpatialDerivatives = Intersection[derList, {"A1", "A2", "A3"}];
    If[firstOrderSpatialDerivatives =!= {}, 
        t = AbsoluteTiming[A1A2A3 = Chop[Transpose[CompiledSHA1A2A3[lMax, firstOrderSpatialDerivatives][Rn, Transpose[sGradientXYZ]], {2, 3, 1}]];][[1]];
        Print["Derivative(s)" <> StringRiffle[firstOrderSpatialDerivatives, {"(", ", ", ")"}] <> " calculated (" <>ToString[t]<>"s)"];
    ];

    secondOrderSpatialDerivatives = Intersection[derList, {"A1A1", "A1A2", "A1A3", "A2A2", "A2A3", "A3A3"}];
    If[secondOrderSpatialDerivatives =!= {}, 
        t = AbsoluteTiming[A11A12A13A22A23A33 = Chop[Transpose[CompiledSHSpatialHessian[lMax, secondOrderSpatialDerivatives][Rn, Transpose[sHessianXYZ]], {2, 3, 1}]];][[1]];
        Print["Derivative(s)" <> StringRiffle[secondOrderSpatialDerivatives, {"(", ", ", ")"}] <> " calculated (" <>ToString[t]<>"s)"];
    ];

    firstOrderAngularDerivatives = Intersection[derList, {"A4", "A5"}];
    If[firstOrderAngularDerivatives =!= {}, 
        t = AbsoluteTiming[A4A5 = Chop[Transpose[CompiledSHA4A5[lMax, firstOrderAngularDerivatives][Rn, sV], {2, 3, 1}]];][[1]];
        Print["Derivative(s)" <> StringRiffle[firstOrderAngularDerivatives, {"(", ", ", ")"}] <> " calculated (" <>ToString[t]<>"s)"];
    ];
    
    secondOrderAngularDerivatives = Intersection[derList, {"A4A4", "A4A5", "A5A5"}];
    If[secondOrderAngularDerivatives =!= {}, 
        t = AbsoluteTiming[A44A45A55 = Chop[Transpose[CompiledSHA44A45A55[lMax, secondOrderAngularDerivatives][Rn, sV], {2, 3, 1}]];][[1]];
        Print["Derivative(s)" <> StringRiffle[secondOrderAngularDerivatives, {"(", ", ", ")"}] <> " calculated (" <>ToString[t]<>"s)"];
    ];

    mixedDerivatives = Intersection[derList, {"A1A4", "A2A4", "A3A4", "A1A5", "A2A5", "A3A5"}];
    If[mixedDerivatives =!= {}, 
        t = AbsoluteTiming[mixedDer = Transpose[CompiledSHMixedDerivativesFromGradientR3[lMax, mixedDerivatives][Rn, Transpose@sGradientXYZ], {2, 3, 1}];][[1]];
        Print["Derivative(s)" <> StringRiffle[mixedDerivatives, {"(", ", ", ")"}] <> " calculated (" <>ToString[t]<>"s)"];
    ];

    Chop[derList /.Join[
                    If[firstOrderSpatialDerivatives =!= {}, MapThread[#1->#2&, {firstOrderSpatialDerivatives, A1A2A3}], {}], 
                    If[secondOrderSpatialDerivatives =!= {}, MapThread[#1->#2&, {secondOrderSpatialDerivatives, A11A12A13A22A23A33}], {}], 
                    If[firstOrderAngularDerivatives =!= {}, MapThread[#1->#2&, {firstOrderAngularDerivatives, A4A5}], {}], 
                    If[secondOrderAngularDerivatives =!= {}, MapThread[#1->#2&, {secondOrderAngularDerivatives, A44A45A55}], {}], 
                    If[mixedDerivatives =!= {}, MapThread[#1->#2&, {mixedDerivatives, mixedDer}], {}]
                ]]

];

LeftInvariantHessianUsingSphericalHarmonics[oso_/;Head[oso]===Obj3DPositionOrientationData, lMax_,σInternal_] := LeftInvariantHessianUsingSphericalHarmonics[oso["Data","Real"], oso[["Wavelets"]][["OrientationList"]], lMax,σInternal]
LeftInvariantHessianUsingSphericalHarmonics[OS_, orientations_, lMax_,σInternal_] := Module[{Rn, A11,A12,A13,A22,A23,A33, A34, A35, A44, A45, A55, A1, A2,A14,A24,A3,A15,A25,A4,A5, A6, hessian, σSpat, σAng,sV, sGradientXYZ, sHessianXYZ,s,dim,R},        
    {σSpat, σAng} = σInternal;
    (*Initialize rotation matrices*)
    Rn = Map[RotateZY[Sequence@@#]&, orientations];

    (*Calculate spherical harmonic coeffiecients*)
    t = AbsoluteTiming[s = PseudoInverseiDSHT[orientations, lMax, OS];][[1]];
    Print["Spherical Harmonic Transform performed (" <>ToString[t]<>"s)"];
    
    (*Calculate first order and second order spatial derivatives*)
    t = AbsoluteTiming[{sV, sGradientXYZ, sHessianXYZ} = HRegularizeAndPrepareSpatialDerivativesXYZ[s, lMax, σSpat, σAng];][[1]];
    Print["Spatial Derivatives Calculated (" <>ToString[t]<>"s)"];
    
    (*Flatten dimensions*)
    dim = Dimensions[OS];
    R = Flatten[ConstantArray[Rn, dim[[1;;3]]], 2];
    sHessianXYZ = Transpose@Flatten[Transpose[sHessianXYZ, {4, 1, 2, 3, 5}], 2];
    sGradientXYZ = Transpose@Flatten[Transpose[sGradientXYZ, {4, 1, 2, 3, 5}], 2];
    sV = Flatten[sV, 2];
    
    (*Calculate derivatives*)        
    {A11,A12,A13,A22,A23,A33, A34, A35, A44, A45, A55, A1, A2,A4,A5,A14,A24,A3,A15,A25} = HCalculateDerivatives[sV, sGradientXYZ, sHessianXYZ, R, lMax, {"A1A1","A1A2","A1A3","A2A2","A2A3","A3A3", "A3A4", "A3A5", "A4A4", "A4A5", "A5A5", "A1", "A2", "A4", "A5","A1A4","A2A4","A3","A1A5","A2A5"}];
    A6 = 0*A1;
    
    (*Build Hessian*)        
    hessian = Transpose[
        Developer`ToPackedArray@{
            {A11,    A12,    A13,    A14,    A15-A3,    A2        },
            {A12,    A22,    A23,    A24+A3,    A25,    -A1        },
            {A13,    A23,    A33,    A34-A2,    A35+A1, A6        },
            {A14,    A24,    A34,    A44,    A45,    A5        },
            {A15,    A25,     A35,     A45,     A55,    -A4        },
            {A6,    A6,        A6,        A6,        A6,        A6        }
        },{3,4,1,2}
    ];
    ArrayReshape[hessian,Join[dim,{6,6}]]
]

HesAToHesB = Compile[{{hes, _Real, 2}, {R, _Real, 2}}, R.hes.Transpose[R], RuntimeAttributes{Listable}];

Structureness[os_, RGaugeFrame_, σInternal_,lMax_] := Module[{hes,hesB,B11,B22,B33,B44,B55,B66,str},
    
    hes = LeftInvariantHessianUsingSphericalHarmonics[os["FullData","Real"], os[["Wavelets"]]["FullOrientationList"], lMax,σInternal];
    hesB = HesAToHesB[hes,RGaugeFrame];
    {B11,B22,B33,B44,B55,B66} = Table[hesB〚All,All,All,All,i,i〛,{i,1,6}];
    str = -(B11+B22+B44+B55+B66);
    str
]

Left Invariant Diffusion in the direction of the Gauge Frames

ExponentialCurveSE3PositionAndOrientationAtgC = Compile[{{c, _Real, 1}, {x, _Real, 1}, {R, _Real, 2}, {t, _Real}}, 
    Module[{q, c1, Ω, x1, n1}, 
        q = Norm[c[[4;;6]]];
        Ω = ({{0, -c[[6]], c[[5]]}, {c[[6]], 0, -c[[4]]}, {-c[[5]], c[[4]], 0}});
        c1 = {c[[1]], c[[2]], c[[3]]};
        x1 = If[q == 0, t*c1, t*c1+(1-Cos[q*t])/q^2 Ω.c1+(t/q^2-Sin[q*t]/q^3) Ω.Ω.c1];
        n1 = If[q == 0, {0, 0, 1}, {0, 0, 1}+Sin[q*t]/q( Ω[[All, 3]])+(1-Cos[q*t])/q^2 (Ω.Ω)[[All, 3]]];
        {R.x1+x, R.n1}
    ], RuntimeAttributes->Listable
];

SpatialInterpolationFactors = Compile[{{pos, _Real, 1}}, 
    Module[{interPolationFactors, rest, x,y,z}, 
        rest = pos - Floor[pos];
        {x, y, z} = rest;
        interPolationFactors = {
            (1-x) (1-y) (1-z), 
            (1-x) (1-y) z, 
            (1-x) y (1-z), 
            (1-x) y z, 
            x (1-y) (1-z), 
            x (1-y) z, 
            x y (1-z), 
            x y z
        }
    ], RuntimeAttributes -> {Listable}
];

Inverse3x3C = Compile[{{mat3x3, _Real, 2}}, 
    Module[{a, b, c, d, e, f, g, h, i, det}, 
        {{a, b, c}, {d, e, f}, {g, h, i}} = mat3x3;
        det = -c e g+b f g+c d h-a f h-b d i+a e i;
        {
            {(-f h+e i)/det, (c h-b i)/det, (-c e+b f)/det}, 
            {(f g-d i)/det, (-c g+a i)/det, (c d-a f)/det}, 
            {(-e g+d h)/det, (b g-a h)/det, (-b d+a e)/det}
        }
    ]
];

Triangles[orientations_List, topology_List] := Triangles[orientations, topology] = Map[orientations[[#]]&, topology, {2}];

JoinedTriangles[orientations_List, topology_List] := JoinedTriangles[orientations, topology] = Position[topology, elem_/;elem == #][[All, 1]]&/@Range[Length[orientations]];

CrossProducts[triangles_List] := CrossProducts[triangles] = Map[{#[[1]]#[[2]], #[[2]]#[[3]], #[[3]]#[[1]]}&, triangles, {1}];

Sides[triangles_List] := Sides[triangles] = Map[{Sign[#[[1]]#[[2]].#[[3]]], Sign[#[[2]]#[[3]].#[[1]]], Sign[#[[3]]#[[1]].#[[2]]]}&, triangles, {1}];

IndexNearestOrientationC = Compile[{{n, _Real, 1}, {orientations, _Real, 2}}, 
    Module[{nearest, i, d, mind}, 
        d = Map[Norm[# - n] &, orientations];
        mind = Min[d];
        i = Position[d, mind][[1, 1]];
        If[i >=  1, i, 1]]
];

AngularInterpolationFactors[or_List, orientations_, topology_] := Module[{triangles = Triangles[orientations, topology], joinedtriangles = JoinedTriangles[orientations, topology], crossproducts, sides, maxLength}, 
    maxLength = Max[Map[Length, joinedtriangles]];
    joinedtriangles = Map[ArrayPad[#, {0, maxLength-Length[#]}]&, joinedtriangles];
    crossproducts = CrossProducts[triangles];
    sides = N[Sides[triangles]];
    AngularInterpolationFactorsPreComputedC[or, orientations, topology, crossproducts, triangles, joinedtriangles, sides]
]

AngularInterpolationFactorsPreComputedC = Compile[{{or, _Real, 1}, {orientations, _Real, 2}, {topology, _Integer, 2}, {crossproducts, _Real, 3}, {triangles, _Real, 3}, {joinedtriangles, _Integer, 2}, {sides, _Real, 2}}, 
    Module[{range, e1, e2, u, row, triangle, triangleid, pos, iNearest, indices, weights}, 
        iNearest = IndexNearestOrientationC[or, orientations];
        If[MemberQ[orientations, or], 
            indices = {iNearest, iNearest, iNearest};
            weights = {1., 0., 0.};
            , 
            (*Select possible triangles by finding triangles connected to the nearest orientation*)
            range = joinedtriangles[[iNearest]];
            
            (*Select triangle*)
            pos = 1;
            While[pos<Length[range]&&range[[pos]] != 0 && (UnitStep[crossproducts[[range[[pos]], All, #]].or]*sides[[range[[pos]], #]]&/@Range[3]) != {1, 1, 1}, 
                pos++];
            triangleid = range[[pos]];
            
            (*Calculate weights and indices*)
            triangle = triangles[[triangleid]];
            e1 = triangle[[2]]-triangle[[1]];
            e2 = triangle[[3]]-triangle[[1]];
            u = Inverse3x3C[{e1, e2, or}].(or-triangle[[1]]);
            weights = {(1.-u[[1]]-u[[2]]), u[[1]], u[[2]]};
            indices = topology[[triangleid]];
        ];
        {indices, weights}
    ], RuntimeAttributes->Listable, CompilationOptions -> {"InlineExternalDefinitions" -> True, "InlineCompiledFunctions" -> True}
];

PrepareInterpolationStep[c_, h_, orientations_, topology_] := Module[{Rn, step, pos, or, floorsS, factS, res, factA, indicesA}, 
    (*Initialize rotation matrices*)
    Rn = Map[RotateZY[Sequence@@#]&, orientations];

    (*Determine orientation and position of step over exponential curve*)
    step = MapIndexed[ExponentialCurveSE3PositionAndOrientationAtgC[#1, #2[[1;;3]], Rn, h]&, c, {3}];
    
    or = step[[All, All, All, All, 2]];
    pos = step[[All, All, All, All, 1]]+1;

    (*Calculate floor and spatial interpolation stencil*)
    floorsS = Floor[pos];
    factS = SpatialInterpolationFactors[pos];

    (*Calculate indices of the vertices used in interpolation and calculate angular interpolation stencil*)
    res = AngularInterpolationFactors[or, orientations, topology];
    factA = res[[All, All, All, All, 2]];
    indicesA = Round[res[[All, All, All, All, 1]]];

    {floorsS, indicesA, factS, factA}
]

PrepareInterpolation[c_, h_, orientations_, topology_] := Module[{intFactorsAndIndicesF, intFactorsAndIndicesB}, 
    intFactorsAndIndicesF = PrepareInterpolationStep[c, h, orientations, topology];
    intFactorsAndIndicesB = PrepareInterpolationStep[c, -h, orientations, topology];
    {intFactorsAndIndicesF, intFactorsAndIndicesB}
]

ExplicitIntegration[U_, Q_, t_List, Δt_] := Module[{output = U, i = 1,tEnd = t〚-1〛,tCurrent,outputAll = {}},
    For[tCurrent = 0, tCurrent ≤ tEnd, tCurrent+= Δt,
        (*See if current endtime has been reached*)
        If[tCurrent + $MachineEpsilon ≥ t〚i〛, AppendTo[outputAll,output]; i++];
        
        (*Perform diffusion step*)
        output = output+Δt*Q[output];
        ];
    AppendTo[outputAll,output];
    Return[outputAll]
]

ExplicitIntegration[U_, Q_, t_, Δt_] := ExplicitIntegration[U, Q, {t}, Δt]〚1〛;

AngularStepsizeCheck[intf_] := Module[{allInterpolationContainIndex}, 
    allInterpolationContainIndex = Map[Intersection[#]&, Transpose[Flatten[Transpose[Flatten[intf[[All, 2]], 3], {2, 3, 1}], 1]]];
    And@@MapIndexed[MemberQ[#1, #2[[1]]]&, allInterpolationContainIndex]
]

LinearInterpolateOSFactorsPrecalculatedC = Compile[{{data, _Real, 4}, {floors, _Integer, 1}, {indices, _Integer, 1}, {interPolationFactors, _Real, 1}, {interPolationFactorsAngular, _Real, 1}}, 
    Module[{dataDim, dataAtIntegerIndices, spatInterPolatedData, interPolatedData, rest, res}, 
        dataAtIntegerIndices = Flatten[data[[floors[[1]] ;; floors[[1]] + 1, floors[[2]] ;; floors[[2]] + 1, floors[[3]] ;; floors[[3]] + 1, indices]], 2];
        spatInterPolatedData = Total[dataAtIntegerIndices*interPolationFactors];
        interPolatedData = spatInterPolatedData.interPolationFactorsAngular;
        interPolatedData
    ], RuntimeAttributes -> {Listable}, CompilationOptions -> {"InlineExternalDefinitions" -> True, "InlineCompiledFunctions" -> True}
];

GroupDifferenceExponentialCurvePrepared[U_, preparedInterpolationFactors_, h_] := Module[{Upad, Uf, Ub}, 
    Upad = ArrayPad[U, 1, "Reversed"][[All, All, All, 2;;-2]];
    Uf = LinearInterpolateOSFactorsPrecalculatedC[Upad, Sequence@@preparedInterpolationFactors[[1]]];
    Ub = LinearInterpolateOSFactorsPrecalculatedC[Upad, Sequence@@preparedInterpolationFactors[[2]]];
    (Ub-2U+Uf)/(2h)^2
]

LIDGAUGE::stabilitycriterion = "Old stability criterion (for non gauge case) is still used";
LIDGAUGE::spatstepLarge = "Spatial stepsize is larger than 1 pixel";
LIDGAUGE::angstepLarge = "Angular stepsize too large (point used in interpolation does not lie in connected mesh triangle)";
Options[LeftInvariantDiffusionGauge] = {DeltaT->Automatic, "StepSize"->0.8};

LeftInvariantDiffusionGauge[os_/;Head[os]===Obj3DPositionOrientationData, RG_, D11_, D33_, D44_, t_/;NumberQ[t], muFrame_, opt___:OptionsPattern[]] := LeftInvariantDiffusionGauge[os, RG, D11, D33, D44, {t}, muFrame, opt]〚1〛;
LeftInvariantDiffusionGauge[os_/;Head[os]===Obj3DPositionOrientationData, RG_, D11_, D33_, D44_, t_List, muFrame_, opt___:OptionsPattern[]] := Module[
    {dataDif, orientations = os[["Wavelets"]]["FullOrientationList"], topology = os[["Wavelets"]]["FullTopology"]},
    dataDif = LeftInvariantDiffusionGauge[os["FullData","Real"], orientations, topology, RG, D11, D33, D44, t, muFrame, opt];
    Affix[os,"Data" #]&/@dataDif
]

LeftInvariantDiffusionGauge[U_List, orientations_List, topology_List, RG_, D11_, D33_, D44_, t_, muFrame_, OptionsPattern[]] := 
Module[{integrator, Q, StabilityCriterion, dt = OptionValue["DeltaT"], spatialstepsizeMax, angularstepsizeMax, h = OptionValue["StepSize"], intF1, intF2, intF3, intF4, intF5, intF6}, 
    (*Message[LIDGAUGE::stabilitycriterion];*)
    If[(h*Max@Sqrt[(Total[(Flatten[RG, Depth[RG]-3][[All, 1;;3]])^2, {2}])])>0.999, Message[LIDGAUGE::spatstepLarge]];
    
    (*mu-orthogonal frame*)
    spatialstepsizeMax = h/muFrame; (*Stepsize in pixels for pure spatial step*)
    angularstepsizeMax = h;(*Angular stepsize for pure angular step*)
    
    (*(*Other frame*)
    spatialstepsizeMax = h*Max@Sqrt[(Total[(Flatten[RG, Depth[RG]-3]〚All, 1;;3〛)^2, {2}])]; 
    angularstepsizeMax = h*Max@Sqrt[(Total[(Flatten[RG, Depth[RG]-3]〚All, 4;;6〛)^2, {2}])];*)
    
    (*Set dt to be equal to the stability bound if it is not specified in advance, otherwise use user input value*)
    StabilityCriterion = 1/((2 Max[D33]/muFrame^2 + 4 Max[D11]/muFrame^2)/spatialstepsizeMax^2+(4Max[D44])/angularstepsizeMax^2);
    If[dt == Automatic, dt = StabilityCriterion];
    
    (*Print grid with stepsizes*)
    Print[
    Grid[{{"Maximum Spatial stepsize", spatialstepsizeMax}, 
        {"Maximum Angular stepsize", angularstepsizeMax}, 
        {"Timestep", dt}}]];
    
    (*Choose integrator*)
    integrator = ExplicitIntegration[#1, #2, #3, #4]&;
    If[dt>StabilityCriterion && integrator === (ExplicitIntegration[#1, #2, #3, #4]&), 
        Message[ContourEnhancement::instable, StabilityCriterion]
    ];
    
    (*Prepare interpolationfactors*)
    intF1 = If[D11 =!= 0, (PrepareInterpolation[RG[[All, All, All, All, 1]], h, orientations, topology]), 0];
    intF2 = If[D11 =!= 0, (PrepareInterpolation[RG[[All, All, All, All, 2]], h, orientations, topology]), 0];
    intF3 = If[D33 =!= 0, (PrepareInterpolation[RG[[All, All, All, All, 3]], h, orientations, topology]), 0];
    intF4 = If[D44 =!= 0, (PrepareInterpolation[RG[[All, All, All, All, 4]], h, orientations, topology]), 0];
    intF5 = If[D44 =!= 0, (PrepareInterpolation[RG[[All, All, All, All, 5]], h, orientations, topology]), 0];
    intF6 = If[D44 =!= 0, (PrepareInterpolation[RG[[All, All, All, All, 6]], h, orientations, topology]), 0];
    (*If[!And@@Map[AngularStepsizeCheck, {intF1, intF2, intF3, intF4, intF5, intF6}], Message[LIDGAUGE::angstepLarge]];*)
    
    (*Form the Q-function needed by the integrators. Note this is only done symbolically*)
    Q[W_] := (
            If[D11  === 0, 0, D11 (GroupDifferenceExponentialCurvePrepared[W, intF1, h]
                                + GroupDifferenceExponentialCurvePrepared[W, intF2, h])
            ]
            + If[D33  === 0, 0, D33 GroupDifferenceExponentialCurvePrepared[W, intF3, h]]
            + If[D44  === 0, 0, D44 (GroupDifferenceExponentialCurvePrepared[W, intF4, h]
                                + GroupDifferenceExponentialCurvePrepared[W, intF5, h]
                                + GroupDifferenceExponentialCurvePrepared[W, intF6, h])
            ]);
            
    (*Do actual diffusion and return result*)
    Return[integrator[U, Q, t, dt]]
]

Left Invariant Gradient

LeftInvariantGradientUsingSphericalHarmonics[oso_/;Head[oso]===Obj3DPositionOrientationData, lMax_,σInternal_] := LeftInvariantGradientUsingSphericalHarmonics[oso["FullData","Real"], oso[["Wavelets"]]["FullOrientationList"], lMax,σInternal];
LeftInvariantGradientUsingSphericalHarmonics[OS_, orientations_, lMax_,σInternal_] := Module[{Rn, A1,A2,A3,A4,A5,A6, hessian, σSpat, σAng,sV, sGradientXYZ, sHessianXYZ,s,dim,R},        
    {σSpat, σAng} = σInternal;
    (*Initialize rotation matrices*)
    Rn = Map[RotateZY[Sequence@@#]&, orientations];

    (*Calculate spherical harmonic coeffiecients*)
    t = AbsoluteTiming[s = PseudoInverseiDSHT[orientations, lMax, OS];][[1]];
    Print["Spherical Harmonic Transform performed (" <>ToString[t]<>"s)"];
    
    (*Calculate first order and second order spatial derivatives*)
    t = AbsoluteTiming[{sV, sGradientXYZ, sHessianXYZ} = HRegularizeAndPrepareSpatialDerivativesXYZ[s, lMax, σSpat, σAng];][[1]];
    Print["Spatial Derivatives Calculated (" <>ToString[t]<>"s)"];
    
    (*Flatten dimensions*)
    dim = Dimensions[OS];
    R = Flatten[ConstantArray[Rn, dim[[1;;3]]], 2];
    sHessianXYZ = Transpose@Flatten[Transpose[sHessianXYZ, {4, 1, 2, 3, 5}], 2];
    sGradientXYZ = Transpose@Flatten[Transpose[sGradientXYZ, {4, 1, 2, 3, 5}], 2];
    sV = Flatten[sV, 2];
    
    (*Calculate derivatives*)        
    {A1,A2,A3,A4,A5} = HCalculateDerivatives[sV, sGradientXYZ, sHessianXYZ, R, lMax, {"A1","A2","A3","A4","A5"}];
    A6 = 0*A1;
    
    (*Build Hessian*)        
    hessian = Transpose[Developer`ToPackedArray@{A1,A2,A3,A4,A5,A6},{3,1,2}];
    ArrayReshape[hessian,Join[dim,{6}]]
]

GradAToGradB = Compile[{{grad, _Real, 1}, {R, _Real, 2}}, R.grad, RuntimeAttributes{Listable}];

V3D Import and Export Functions (for AIXIA viewer)

LoadAndPartitionDataV3D[if_, zdim_,ydim_,xdim_,bitDepth_] := Module[{data, type = Switch[bitDepth, 16, "UnsignedInteger16", 32,"UnsignedInteger32",_,"UnsignedInteger16"]},
    data = ArrayReshape[
        BinaryReadList[if, type, zdim*ydim*xdim], 
        {zdim, ydim, xdim}];
    data
]

LoadDataV3D[fn_] := Module[{if,isFiltered,zdim, ydim, xdim,data, voxelSize, bitDepth,percentage, centerPosition, metaData, version},
    if = OpenRead[fn, BinaryFormat -> True];
    
    (*Import version*)
    version = StringJoin@BinaryReadList[if, "Character8", 40];
    
    (*Import data dimension*)
    {xdim, ydim, zdim} = BinaryReadList[if, "UnsignedInteger32", 3];
    
    (*Import voxels size*)
    voxelSize = BinaryReadList[if, "Real64", 3];
    
    (*Import data center*)
    centerPosition = BinaryReadList[if, "Real64", 3];

    (*Import BitDepth*)
    bitDepth = BinaryRead[if, "UnsignedInteger32"];

    (*Import percentage ???? Nobody knows what it is ???? *)
    percentage = BinaryRead[if, "UnsignedInteger32"];
    
    (*Import slope and intercept*)
    isFiltered = BinaryRead[if, "UnsignedInteger32"];
    
    (*Import data and partition*)
    data = LoadAndPartitionDataV3D[if, zdim, ydim, xdim,bitDepth];
    
    (*Close file and return results*)
    Close[if];
    
    metaData = {"Version"-> version, "BitDepth" -> bitDepth, "VoxelSize"-> voxelSize, "Center"-> centerPosition, "IsFiltered"-> isFiltered, "Percentage" -> percentage, "Dimensions" -> {zdim, ydim, xdim}};
    {data,metaData}
]

SaveDataV3D::mtd = "Metadata does not contain all required fields or no meta data is given. Default values are chosen for the following fields: `1`";
SaveDataV3D[fn_, data_, metaData_:{}] := Module[{if,isFiltered, dataDimensions, voxelSize, bitDepth,percentage, centerPosition, version, type, allKeys, givenFields, missingFields,defaultValues},
    (*Open write file*)
    if = OpenWrite[fn, BinaryFormat -> True];
    
    (*Check which keys are present in metadata*)
    allKeys = {"Version", "BitDepth","VoxelSize","Center","IsFiltered","Percentage"};
    givenFields = Intersection[allKeys,Keys@metaData];
    missingFields = Complement[allKeys,Keys@metaData];
    defaultValues = {
        "Version" -> StringJoin["3D-RA R6.1",StringRepeat[" ", 40 - StringLength["3D-RA R6.1"]]],
        "BitDepth" -> 0,
        "VoxelSize" -> {1.,1.,1.},
        "Center" -> {0,0,0},
        "IsFiltered" -> 0,
        "Percentage" -> 100
    };
    
    If[missingFields =!= {}, Message[SaveDataV3D::mtd,missingFields]];
    
    {version, bitDepth,voxelSize,centerPosition,isFiltered,percentage} = {"Version", "BitDepth","VoxelSize","Center","IsFiltered","Percentage"}/.metaData/.defaultValues;
    dataDimensions = Reverse@Dimensions[data];
    
    (*Save version*)
    BinaryWrite[if, version, "Character8"];
    
    (*Save data dimension*)
    BinaryWrite[if, dataDimensions, "UnsignedInteger32"];
    
    (*Save voxels size*)
    BinaryWrite[if, voxelSize, "Real64"];
    
    (*Save data center*)
    BinaryWrite[if, centerPosition, "Real64"];

    (*Save BitDepth*)
    BinaryWrite[if, bitDepth, "UnsignedInteger32"];

    (*Save percentage ???? Nobody knows what it is ???? *)
    BinaryWrite[if, percentage, "UnsignedInteger32"];
    
    (*Save slope and intercept*)
    BinaryWrite[if, isFiltered, "UnsignedInteger32"];
    
    (*Save data*)
    type = Switch[bitDepth, 16, "UnsignedInteger16", 32,"UnsignedInteger32",_,"UnsignedInteger16"];
    BinaryWrite[if, data, type];
    
    (*Close file and return results*)
    Close[if];
]