function figvec = satlib_plot_ungrowth_data(varargin)
% function satlib_plot_ungrowth_data(varargin)
%
% Shows the areal growth real and anisotropy for ungrowth data
%
% example 1
% show_region_areal_growth;
%
% example 2
% show_region_areal_growth('ungrowthfilename', 'Staged\Data\UngrowthData\Ungrowth_aih_001.mat', 'scale', 0);
%
% Dr. A. I. Hanna (2007)
data = parseargs(varargin);
ungrown_grids = data.ungrown_grids;

boundary = tri2boundary(ungrown_grids{1}.tri_info.T, ungrown_grids{1}.tri_info.pts);



ug = data.ungrown_grids;
N = length(ug);
lims = [realmax, realmin];
values = {};
for i=2:N
    if ~isfield(ug{i}, 'deltaSt')
        fprintf('WARNING: could not find a deltaSt for this stage, setting to 20 for default\n');
        T = 20;
    else
        T = ug{i}.deltaSt;
    end
    fprintf('Using %3.3f as the time in hours between stages %d and %d\n', T, i-1, i);
    if strcmpi(data.gptype, 'mesh')
        growthRates = growthParameters2GrowthRates(satio_meshes2GrowthData(ug{i-1}, ug{i}, T), T);
    else
        growthRates = growthParameters2GrowthRates(ug{i}.ungrowth_tensors, T);
    end
    [vals, str] = getSchema(growthRates, data.method);
    if min(vals)<lims(1)
        lims(1) = min(vals);
    end
    if max(vals)>lims(2)
        lims(2) = max(vals);
    end
    values = cat(1, values, vals);
end
[map, lims] = getColormap(data.colormap, data.flipcolmap, lims);

if ~isempty(data.userlims)
    lims = data.userlims;
end
colours = cell(length(values)+1, 1);
if ~strcmpi(data.method, 'uniform')
    for i=1:length(values)
        vals = values{i};
        vals = (vals - lims(1))./(lims(2)-lims(1));
        ind = vals<0;
        vals(ind) = 0;
        ind = vals>1;
        vals(ind) = 1;
        
        cols = map(floor(254*vals) + 1,:);
        if strcmpi(data.colormap, 'lab')
            cols(ind,1) = 1;
            cols(ind,2) = 1;
            cols(ind,3) = 1;
        end
        colours{i+1} = cols;
    end
    colours{1} =ones(size(ug{1}.ungrowth_tensors,3), 3);
    data.facecolor = colours;
end
figvec = satlib_drawRegions('grids', ug, 'colours', data.facecolor, 'scale', data.scale, ...
    'title', str, 'EdgeColor', data.edgecolor, 'ShowTensors', data.showtensors, ...
    'outlinecolour', data.outlinecol, 'scale', data.scale, 'tensorcol', data.tensorcol, ...
    'colormap', map, 'collims', lims, 'linewidth', data.linewidth, 'boundary', boundary, ...
    'tensorscalefactor', data.tensorscalefactor, 'tensorthreshold', data.tensorthreshold, 'tensortype', data.tensortype);
return;
%%
function range = extendToZero( range )
    if length(range)==1
        if range <= 0
            range = [ range 0];
        else
            range = [0 range];
        end
    else
        if range(1) >= 0
            range(1) = 0;
        end
        if range(2) <= 0
            range(2) = 0;
        end
    end
return
%%
function c = posnegMap( range, zerowhite, negpart, pospart )
    if nargin==0, range = [-1 1]; end
    range = extendToZero( range );
    if range(1)==range(2)
        c = [1 1 1
             1 1 1];
        return;
    end
    
    if nargin < 2
        zerowhite = true;
    end
    nsteps = 254;
    if zerowhite
        negfrac = range(1)/(range(1)-range(2));
        posfrac = range(2)/(range(2)-range(1));
        nneg = ceil(negfrac*nsteps);
        npos = ceil(posfrac*nsteps);
        maxfrac = max(negfrac,posfrac);
        negfrac = negfrac/maxfrac;
        posfrac = posfrac/maxfrac;
        if nneg==0
            posmap = makeCmap( pospart, npos+1, posfrac );
            c = posmap;
        elseif npos==0
            negmap = makeCmap( negpart, nneg, negfrac );
            c = negmap(size(negmap,1):-1:1,:);
        else
            negmap = makeCmap( negpart, nneg, negfrac );
            posmap = makeCmap( pospart, npos, posfrac );
            c = [ negmap(size(negmap,1):-1:1,:);
                  posmap(2:size(posmap,1),:) ];
        end
    else
        negpos = [ negpart( size(negpart,1):-1:1, : ); pospart ]; 
        c = makeCmap( negpos, nsteps, 1 );
    end
return
%%
function cmap = makeCmap( colors, nsteps, frac )
    if nsteps <= 0
        cmap = zeros(0,3);
        return;
    end
    cmap = zeros( nsteps, 3 );
    ncsteps = size(colors,1)-1;
    a = (0:(nsteps))*frac;
    b = a*(ncsteps/nsteps);
    c = floor(b);
    d = b-c;
    e = 1-d;
    for i=1:nsteps
        ci = c(i);
        cmap(i,:) = colors(ci+1,:)*e(i) + colors(ci+2,:)*d(i);
    end
return

%%
function [c,range] = positiveRainbow( range )
    negpart = [ 1 1 1       % white
                0.9 0 1       % purple
                0 0 0 ];    % black

    pospart = [ 1 1 1       % white
                0 0 1       % blue
                0 0.7 1     % bluish cyan
                0 1 1       % cyan
                0 1 0.7     % greenish cyan
                0 1 0       % green
                0.75 1 0    % chartreuse
                1 1 0       % yellow
                1 0.875 0   % yellow-orange
                1 0.5 0     % orange
                1 0 0       % red
                0.6 0 0 ];  % darker red

    range = extendToZero( range );
    c = posnegMap( range, true, negpart, pospart );
return

%%
function [map, lims] = getColormap(cmapname, flipcolmap, lims)
if strcmpi(cmapname, 'lab')
    map = jet(255);
    if lims(1)<0
        lims(1) = 0;
    end
elseif strcmpi(cmapname, 'lab2')
    [map, range] = positiveRainbow(lims);
    if lims(1)<0
        lims(1) = 0;
    end
else
    
    map = feval(cmapname, 255);
end
if  flipcolmap
    map = flipud(map);
end
return
%%
function [vals, str] = getSchema(growthRates, method)
vals = 1;
str = 'Uniform';
switch lower(method)
    case 'smaj'
        vals = [growthRates.Smaj];
        str = 'Smaj';
    case 'smin'
        vals = [growthRates.Smin];
        str = 'Smin';
    case 'sarea'
        vals = [growthRates.Sarea];
        str = 'Sarea';
    case 'kmaj'
        vals = [growthRates.Kmaj];
        str = 'Kmaj';
    case 'kmin'
        vals = [growthRates.Kmin];
        str = 'Kmin';
    case 'karea'
        vals = [growthRates.Karea];
        str = 'Karea';
    case 'kmaj/kmin'
          vals = [growthRates.Kanis];
        str = 'Kmaj/Kmin';
    case 'rmaj'
        vals = [growthRates.Rmaj];
        str = 'Rmaj';
    case 'rmin'
        vals = [growthRates.Rmin];
        str = 'Rmin';
    case 'rarea'
        vals = [growthRates.Rarea];
        str = 'Rarea';
    case 'dmaj'
        vals = [growthRates.Dmaj];
        str = 'Dmaj';
    case 'dmin'
        vals = [growthRates.Dmin];
        str = 'Dmin';
    case 'darea'
        vals = [growthRates.Darea];
        str = 'Darea';
end
return
%%
%
function data = parseargs(varargin)
varargin = varargin{1};
if mod(length(varargin),2) ~= 0
    % input args have not com in pairs, woe is me
    error('Arguments to show_region_areal_growth must come param/value in pairs.')
end
data.gptype = 'sector';
data.tensortype = [1 1];
data.ungrowthfilename = [];
data.scale = 1;
data.linewidth = 1;
data.outlinecol = 'k';
data.showmesh = 1;
data.showtensors = 1;
data.edgecolor = 'k';
data.facecolor = [0 1 0];
data.tensorcol = [0 0 0];
data.method = 'uniform';
data.colormap = 'hsv';
data.userlims = [];
data.ungrown_grids = {};
data.tensorscalefactor = 1;
data.tensorthreshold = 0;
data.tensoruniscale = 0;
data.flipcolmap = 0;
for i=1:2:length(varargin)
    switch lower(varargin{i})
        case 'ungrown_grids'
            data.ungrown_grids = varargin{i+1};
        case 'gptype'
            data.gptype = varargin{i+1};
        case 'tensortype'
            data.tensortype = varargin{i+1};
        case 'tensoruniscale'
            data.tensoruniscale = varargin{i+1};
        case 'tensorscalefactor'
            data.tensorscalefactor = varargin{i+1};
        case 'tensorthreshold'
            data.tensorthreshold = varargin{i+1};
        case 'ungrowthfilename'
            data.ungrowthfilename = varargin{i+1};
        case 'scale'
            data.scale = varargin{i+1};
        case 'outlinecolour'
            data.outlinecol = varargin{i+1};
        case 'showmesh'
            data.showmesh = varargin{i+1};
        case 'showtensors'
            data.showtensors = varargin{i+1};
        case 'edgecolor'
            data.edgecolor = varargin{i+1};
        case 'facecolor'
            data.facecolor = varargin{i+1};
        case 'tensorcol'
            data.tensorcol = varargin{i+1};
        case 'method'
            data.method = varargin{i+1};
        case 'colormap'
            data.colormap = varargin{i+1};
        case 'flipcolormap'
            data.flipcolmap = varargin{i+1};
        case 'userlimits'
            data.userlims = varargin{i+1};
        case 'linewidth'
            data.linewidth = varargin{i+1};
        otherwise
            error(['Unknown parameter name passed to ', mfilename, ' show_region_areal_growth.  Name was ' varargin{i}])
    end
end

if isempty(data.ungrown_grids)
    if isempty(data.ungrowthfilename)
        [filename, pathname] = uigetfile('*.mat', 'Pick an Ungrowth File');
        if isequal(filename,0) || isequal(pathname,0)
            disp('User pressed cancel')
            return;
        end
        data.ungrowthfilename = [pathname, filesep, filename];
    end
end
