close all;
clear all;
params = get_params();
startAngles ={275,310,325,355};

penetrationSpan = 2:600;

dataFolderBase = '../data/';
dataSet = 'SAS3';

dataFolder = fullfile(dataFolderBase,dataSet);
resultsFolder = fullfile('../results/',dataSet);

fprintf('Running dataSet %s\n',dataSet);
fieldfilename = 'field2d.mat';
fieldFilePath = fullfile(dataFolder,fieldfilename);
if ~exist(fieldFilePath,'file')
	fprintf('%s file not found\n',fieldfilename);
	fprintf('Please run "generateOrientationField" script for datset %s to generate it\n',dataSet);
end
data2d = load(fieldFilePath);
field = data2d.field2D;
mask = data2d.mask2D;
background = mat2gray(data2d.wga2D);

[helixAngles,penMask,X] = fieldToHelixAngle2(field,mask);

plotsFolder = fullfile(resultsFolder,'/sectorPlots/');
if not(isfolder(plotsFolder))
	mkdir(plotsFolder);
    fprintf('Created Folder %s\n',plotsFolder)
else
    fprintf('Folder %s already exists.. old plots may be overwritten\n',plotsFolder);
end

for k = 1:length(startAngles) 
	fprintf("Plotting angle (%d/%d)\n",k,length(startAngles));
	sectorStartAngles = startAngles{k}; %degrees
	startAngleString = [];
	angularWidth = 10; % degrees
	sector = [];
	selectedIndices = [];
    for s = 1:length(sectorStartAngles)
        startAngle = sectorStartAngles(s);
        [currentSector, currentSelectedIndices] = extractHelixAngleSector(helixAngles,startAngle,angularWidth);
        currentSector = rad2deg(medianGuided_angle_rectify(currentSector));
        sector = [sector;currentSector];
        selectedIndices = [selectedIndices, currentSelectedIndices];
        startAngleString = [startAngleString,num2str(startAngle,'%03d'),'-'];
    end

    clear startAngle;
    sectorMask = false(size(penMask));
    s = 1;
    tmphelixAngles = nan*ones(size(helixAngles));
    for r = selectedIndices
       sectorMask(penMask == r) = true;
       if params.useColormappedWedge
            tmphelixAngles(r,:) = deg2rad(sector(s,:));
            s = s+1;
       end
    end
    sectorMask = imdilate(sectorMask,strel('disk',5)) & mask;
    sector = sector(:,penetrationSpan);
    xAxis = linspace(0,length(penetrationSpan)*params.voxelSize,length(penetrationSpan));
    fprintf("Extracted sector(s) : %s\n",startAngleString);
    switch(params.averaging_method)
        case 'mean'
            averageSectorPlot = mean(sector,'omitnan');
        case 'median'
            averageSectorPlot = median(sector,'omitnan');
        otherwise
            averageSectorPlot = median(sector,'omitnan');
    end
    rgb = zeros([size(background),3]);
    if params.useColormappedWedge
        [~, cmapcolored] = helixAnglesToHelixMap(tmphelixAngles,X,mask,sectorMask);
    end
    for channel = 1:3
        if params.useColormappedWedge
            rgb(:,:,channel) = imcomplement(background).*(~sectorMask) + cmapcolored(:,:,channel).*(sectorMask);
        else
            rgb(:,:,channel) = imcomplement(background).*(~sectorMask) + background.*(sectorMask).*params.wedgeColor(channel);
        end
    end
    sectorName  = ['sector-' , startAngleString,'w-' , num2str(angularWidth,'%03d')];
%     if params.loadSmooth
%         sectorName = ['smooth-',sectorName];
%     end
    
    if params.useSinglePlotColor
        plotColors = repmat(params.singleColor,[size(sector,1),1]);
    else
        plotColors = parula(size(sector,1));
    end
    fig = figure('Visible',params.display,'Position',[0,0,params.W,params.H],'DefaultAxesFontSize',18);
    ax = subplot(1,1,1);
    xlabel('Penetration Depth (voxel units) ->');
    ylabel('Angle (degrees) -> ');
    
    ylim(ax,[-200,200]);
    hold(ax,'on')
    for r = 1:size(sector,1)
        index = mod(r,size(plotColors,1)) + 1;
        color = [plotColors(index,:),params.plotAlpha];
        plot(ax,xAxis, sector(r,:),'Color',color,'LineWidth',0.75);
    end
    
    plot(ax,xAxis, averageSectorPlot,'Color',params.medianColor,'LineStyle',':','LineWidth',2.75);
    text(ax,params.xzoom+params.w+0.05,0.96,[startAngleString,'^\circ(',num2str(angularWidth),'^\circ)'],...
        'Color','red','FontSize',18,'Units','normalized')
    if params.savePlots
        plotFileName = fullfile(plotsFolder, ['plot-',sectorName,'.eps']);
        exportgraphics(fig,plotFileName,'Resolution',1200);
    end
    subax = axes(fig,'position',[params.xstart params.ystart params.w params.h]);
    imshow(rgb,'Parent',subax);
    
    zoomax = axes(fig,'position',[params.xzoom params.yzoom params.w-0.275 params.h-0.05 ]);
    plot(zoomax,xAxis(params.zoomRegion),averageSectorPlot(params.zoomRegion),'Color',params.medianColor,'LineWidth',2);
    %box(zoomax,'off');
    zoomax.YAxisLocation = 'left';
    zoomax.YColor = 'blue';
    zoomax.XColor = 'blue';
    zoomax.XAxisLocation = 'bottom';
    
    combinedPlotName = ['combined-',sectorName,'.png'];
    if params.savePlots
        combinedPlotFileName = fullfile(plotsFolder,combinedPlotName);
        exportgraphics(fig,combinedPlotFileName,'Resolution',600);
    end
end

function params = get_params()
    params.averaging_method = 'median';%'mean';
    params.W = 1048;
    params.H = 1048;
    
    
    params.w = 0.36;
    params.h = params.w;
    
    params.xstart = 0.905 - params.w;
    params.ystart = 0.115;
    
    params.xzoom=0.2;
    params.yzoom= 0.975 - params.h;
    params.zoomRegion = 1:100;
    
    params.display = 'on';
%    if isMATLABReleaseOlderThan("R2020a","release")
    rel = version('-release');
    cmp = le('2020a',rel);
    if ~all(cmp)
        warning("This version of matlab does not support saving plots")
        params.savePlots = false;
    else
        params.savePlots = true;
    end
    params.medianColor = [1,0,1];
    params.wedgeColor = [1,0,1];
    params.useColormappedWedge =true;%false;
    
    params.voxelSize = 1.98;
    params.useSinglePlotColor = true;
    
%     params.loadSmooth =false;
%     params.smoothFile = 'field2_step_0025.mat';
    
    params.singleColor = [0.8, 0.8, 0.8];
    params.plotAlpha = 0.2;
end