function figvec = satlib_warpRegularMesh(varargin)
% function figvec = satlib_warpRegularMesh(varargin)
%
% Generates and warps a regular grid according to the deformation of the
% supplied ungrown grids.
%
% Uses the motion coherence library to determine the warp. It has been
% proven that for sparse data sets the chosen regularizer in mctlib
% satisfies the two constraints for regularizers put forward by Yuille.
%
% Dr. A. I. Hanna, CMP, 2008
[grids, dx, dy, facecolor, edgecolor, scale, direction, smoothingterm, ...
    showvectors, linewidth, rotval, outlinecol, padval, showclones, clonescale] = parseargs(varargin);

if isempty(grids)
    return
end

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

sb = sortboundary(boundary);


if strcmpi(direction, 'backward')
    grids = flipud(grids);
end
figvec = [];

ngrids = length(grids);
[x, y, xlim, ylim] = getRegularGrid(grids{1}, dx, dy, rotval, grids{1}.tri_info.pts(sb,:), padval);
zscale = min(xlim(2) -xlim(1), ylim(2)-ylim(1));
figh = figure;
axesh = gca;
figvec = cat(1, figvec, figh);

hold on;
sh = surf(x, y, zeros(size(x)), 'FaceAlpha', 1, 'FaceColor', facecolor, 'EdgeColor', edgecolor, 'LineWidth', linewidth);
pts0 = grids{1}.tri_info.pts;
if ~strcmpi(outlinecol, 'none')
    triplotedges(boundary, pts0, linewidth, outlinecol);
end

pts0(sb,:) = [];
theta = linspace(0, 2*pi, 20);
xc = clonescale*cos(theta);
yc = clonescale*sin(theta);
c = [xc(:), yc(:)];
C = [];
for i=1:size(pts0,1)
    rc = c + ones(length(theta),1)*pts0(i,:);
    C = cat(3, C, rc);
    if showclones
        fill(rc(:,1), rc(:,2), 'r');
    end
end

axis image ij;
xlim = get(gca, 'Xlim');
ylim = get(gca, 'Ylim');
for i=2:ngrids
    pts1 = grids{i-1}.tri_info.pts;
    pts2 = grids{i}.tri_info.pts;
    wpts = mctlib_warpptsxy('XX', [x(:), y(:)], 'x0', pts1, 'y0', pts2, 'sigma', 1, 'lambda', smoothingterm);
    x = reshape(wpts(:,1), size(x));
    y = reshape(wpts(:,2), size(y));
    figh = figure;
    axesh = cat(1, axesh, gca);
    hold(gca, 'on');
    figvec = cat(1, figvec, figh);
    sh = surf(x, y, zeros(size(x)), 'FaceAlpha', 1, 'FaceColor', facecolor, 'EdgeColor', edgecolor, 'LineWidth', linewidth);
    if ~strcmpi(outlinecol, 'none')
        triplotedges(boundary, grids{i}.tri_info.pts, linewidth, outlinecol);
    end
    if showvectors
        hold on;
        smoothingterm
        wpts1 = mctlib_warpptsxy('XX', pts1, 'x0', pts1, 'y0', pts2, 'sigma', 1, 'lambda', smoothingterm);
        plot(pts1(:,1), pts1(:,2),'o', 'MarkerSize', 4*linewidth, 'LineWidth', linewidth);
        v = wpts1 - pts1;
        plot(pts2(:,1), pts2(:,2),'go', 'MarkerSize', 4*linewidth, 'LineWidth', linewidth);
        %v = pts2 - pts1;
        quiver(pts1(:,1), pts1(:,2), v(:,1), v(:,2), 0, 'LineWidth', linewidth);
    end
    if showclones
        for c = 1:size(C,3)
            C(:,:,c) = mctlib_warpptsxy('XX', C(:,:,c), 'x0', pts1, 'y0', pts2, 'sigma', 1, 'lambda', smoothingterm);
            fill(C(:,1,c), C(:,2,c), 'r');
        end

    end
    view(2);
    axis image ij;
    [xlim, ylim] = updateLims(xlim, ylim, get(gca, 'Xlim'), get(gca, 'Ylim'))
    drawnow;
end
if scale
    set(axesh, 'XLim', xlim, 'YLim', ylim);
end
return
%%
function [xlim, ylim] = updateLims(xlim, ylim, cxlim, cylim)
if cxlim(1)<xlim(1)
    xlim(1) = cxlim(1);
end
if cxlim(2)>xlim(2)
    xlim(2) = cxlim(2);
end
if cylim(1)<ylim(1)
    ylim(1) = cylim(1);
end
if cylim(2)>ylim(2)
    ylim(2) = cylim(2);
end
return

%%
function sb = sortboundary(b)
sb = b(1,1:2)';
val = b(1,2);
while(val~=sb(1))
    in = b(:,1)==val;
    val = b(in,2);
    sb = cat(1, sb, val);
end
return
%%
%
function [x, y, xlim, ylim] = getRegularGrid(grid, Nx, Ny, rotval, S, padval)
T = grid.tri_info.T;
pts = grid.tri_info.pts;

R = [cosd(rotval), -sind(rotval); sind(rotval), cosd(rotval)];

mu = mean(pts);
pts = pts - ones(size(pts,1),1)*mu;

pts = (R'*pts')';
pts = pts + ones(size(pts,1),1)*mu;

xlim = [min(pts(:,1)), max(pts(:,1))];
m = mean(xlim);
xlim = 2*(xlim-m) + m;
ylim = [min(pts(:,2)), max(pts(:,2))];
m = mean(ylim);
ylim = 2*(ylim-m) + m;
xind = linspace(xlim(1), xlim(2), 2*Nx);
yind = linspace(ylim(1), ylim(2), 2*Ny);
[x, y] = meshgrid(xind, yind);

sz = size(x);
X = [x(:), y(:)];
mu = mean(X);
X = X - ones(size(X,1),1)*mu;
X = (R*X')';
X = X + ones(size(X,1),1)*mu;

x = reshape(X(:,1), sz);
y = reshape(X(:,2), sz);

mu = mean(S);
S = S - ones(size(S,1),1)*mu;
S = padval*S;
S = S + ones(size(S,1),1)*mu;

% only keep the points that are inside the shape
in = ~inpoly([x(:), y(:)], S);
in = reshape(in, size(x));
x(in) = nan;
y(in) = nan;
return
%%
%
function grids = loadGrids(filename)
grids = {};
if ~exist(filename, 'file')
    return;
end
grids = load(filename);
if ~isfield(grids, 'ungrown_grids')
    return;
end
grids = grids.ungrown_grids;
return;
%%
%
function [grids, dx, dy, facecolor, edgecolor, scale, direction, smoothingterm, ...
    showvectors, linewidth, rotval, outlinecol, padval, showclones,clonescale] = parseargs(varargin)
varargin = varargin{1};
if mod(length(varargin),2) ~= 0
    % input args have not com in pairs, woe is me
    error(['Arguments to', mfilename, ' must come param/value in pairs.'])
end
clonescale = 1;
padval = 1;
dx = 20;
dy = 20;
rotval = 0;
grids = {};
facecolor = 'g';
edgecolor = 'k';
showclones =0;
scale = 0;
linewidth = 1;
outlinecol = 'none';
direction = 'forward';
smoothingterm = 0;
showvectors = 0;
for i=1:2:length(varargin)
    switch lower(varargin{i})
        case 'clonescale'
            clonescale = varargin{i+1};
        case 'showclones'
            showclones =varargin{i+1};
        case 'padval'
            padval = varargin{i+1};
        case 'outlinecol'
            outlinecol = varargin{i+1};
        case 'rotval'
            rotval = varargin{i+1};
        case 'linewidth'
            linewidth = varargin{i+1};
        case 'grids'
            grids = varargin{i+1};
        case 'filename'
            filename = varargin{i+1};
        case 'dx'
            dx = varargin{i+1};
        case 'dy'
            dy = varargin{i+1};
        case 'facecolor'
            facecolor = varargin{i+1};
        case 'edgecolor'
            edgecolor = varargin{i+1};
        case 'scale'
            scale = varargin{i+1};
        case 'direction'
            direction = varargin{i+1};
        case 'smoothingterm'
            smoothingterm = varargin{i+1};
        case 'showvectors'
            showvectors = varargin{i+1};
        otherwise
            error(['Unknown parameter name passed to ', mfilename, ' show_region_areal_growth.  Name was ' varargin{i}])
    end
end
if isempty(grids)
    if isempty(filename)
        error('Parameter GRIDS or FILENAME must be supplied');
    else
        grids = loadGrids(filename);
    end
end
return