function [field, mask] = orientationEstimate(wga,rho,sigma)
    if nargin < 2
        rho = 3*[1,1,1]; 
    end
    if nargin < 3
        sigma = 0.5*[1,1,1];
    end
    wga = imgaussfilt3(wga,sigma);
    [U,V,W,~] = structureTensorOrientation(wga,rho);
    mask = imgaussfilt3(wga,1);
    mask = mask > 500;
    field(:,:,:,1) = U;
    field(:,:,:,2) = V;
    field(:,:,:,3) = W;
end

function [U,V,W,m]=structureTensorOrientation(input,rho)
    if length(rho) == 1
        rho = [rho,rho,rho];
    end
    
    firstInput = single(im2double(input));
    fprintf('Computing Gradients\n');
    gradStart = tic;
    [fx, fy, fz] = gradient(firstInput);
    
    Jfxx = imgaussfilt3(fx.^2,rho);
    Jfxy = imgaussfilt3(fx.*fy,rho);
    Jfxz = imgaussfilt3(fx.*fz,rho);
    clear fx;
    Jfyy = imgaussfilt3(fy.*fy,rho);
    Jfyz = imgaussfilt3(fy.*fz,rho);
    clear fy;
    Jfzz = imgaussfilt3(fz.*fz,rho);
    clear fz;
    fprintf('Computing Gradients took %f seconds\n',toc(gradStart));
    
    fprintf('Starting Eigen Vector Decomposition\n This may take a while..\n');
    eigStart = tic;
    [~,~,m,~,~,~,~,~,~,U,V,W] = EigenVectors3D(Jfxx,Jfxy,Jfxz,Jfyy,Jfyz,Jfzz);
                                                    
    fprintf('Eigen Vector Decomposition took %f seconds\n',toc(eigStart));
end