dataRoot = '../data/SAS3/';
fprintf('Generating field for %s\n',dataRoot);

sigma = 4;
if ~exist([dataRoot,'field_sigma_',num2str(sigma),'.mat'],'file')
    fprintf('Field file not found\nStarting field Estimation\n')
    denoised = bfReadVolume(fullfile(dataRoot,'denoised.tif'));
    startTime =tic;
    [field,~] = orientationEstimate(denoised);
    clear denoised;
	[wga,mask3d] = loadrawdata(dataRoot);
    field = field.*cat(4,mask3d,mask3d,mask3d);
    fprintf('running smoothing\nCould take a while longer..Check back after a couple of hours\n')
    field = projectionSmooth(single(field),single(mask3d),sigma);
    save(fullfile(dataRoot,['field_sigma_',num2str(sigma),'.mat']),'field','-v7.3');
    fprintf('Computed orientation field in %.2f seconds\n',toc(startTime));
else
    fprintf('Loading precomputed field\n');
    field = load(fullfile(dataRoot,['field_sigma_',num2str(sigma),'.mat']));
    field = field.field;
    fprintf('Loaded precomputed orientation field for %s\n',dataRoot);
end
if ~exist('wga','var')
	[wga,mask3d] = loadrawdata(dataRoot);
end
mask2D = max(mask3d,[],3);
wga2D = mat2gray(max(wga,[],3));
field2D = maskedCollapseZSmoothing(field,mask3d);
field2D(~isfinite(field2D)) = 0;
save([dataRoot,'field2d.mat'],'field2D','wga2D','mask2D','-v7.3');

function [wga,mask3d] = loadrawdata(dataRoot)
    wga = bfReadVolume(fullfile(dataRoot,'wga.tif'));
    mask3d = imgaussfilt3(single( wga > 20), [30,30,1]) > 0.5;
    mask3d = bwareaopen(mask3d,400);
end