close all;
clear all;
clc;
load('V_CELL');load('V_DATA');load('V_ORIG');

%%
close all; clc;

i = 15;
L_tmp = length(V_cell{i});
Ax_tmp = reshape(V_cell{i}(1, :), L_tmp/6, 6);
Ay_tmp = reshape(V_cell{i}(2, :), L_tmp/6, 6);

Ax = Ax_tmp(:, 1); Ay = Ay_tmp(:, 1);
Bx = Ax_tmp(:, 6); By = Ay_tmp(:, 6);

%Ax = [0:4]; Ay = zeros(1, 5);
%Bx = Ax; By = Ay - 1;

figure, plot(Ax, Ay, 'r-*'), hold on; plot(Bx, By, 'b-*');
for i=1:length(Ax), plot([Ax(i) Bx(i)], [Ay(i) By(i)], 'k', 'linewidth', 1);end


V1x = Ax(1:(end-1)) - Ax(2:end);
V1y = Ay(1:(end-1)) - Ay(2:end);

V2x = Ax(1:(end-1)) - Bx(1:(end-1));
V2y = Ay(1:(end-1)) - By(1:(end-1));

V3x = Ax(1:(end-1)) - Bx(2:end);
V3y = Ay(1:(end-1)) - By(2:end);

Volumes_before = abs(V1x.*V3y - V3x.*V1y)/2 + abs(V2x.*V3y - V3x.*V2y)/2;
Lens_before = sqrt((Ax - Bx).^2 + (Ay - By).^2);
Vo = mean(Volumes); Lo = mean(Lens);
N = 10;
KV = 200; KL = 0; N_steps = 1000;
disp(['KV = ', num2str(KV), ' KL = ', num2str(KL)]);
Vo = Vo;
Lo = Lo;

[Gx, Gy] = CALC_GRADIENT_FULL(Ax, Ay, Bx, By, KV, Vo, KL, Lo);
quiver(Bx, By, Gx, Gy);

[Bx_new, By_new] = E_MIN_EVOLVE(Ax, Ay, Bx, By, KV, Vo, KL, Lo, N_steps);
disp(['E_old = ', num2str(CALC_ENERGY(Ax, Ay, Bx, By, KV, Vo, KL, Lo)), ' E_new = ', num2str(CALC_ENERGY(Ax, Ay, Bx_new, By_new, KV, Vo, KL, Lo))]);

plot(Bx_new, By_new, 'k-*');
for i=1:length(Ax), plot([Ax(i) Bx_new(i)], [Ay(i) By_new(i)], 'k', 'linewidth', 1);end

V1x = Ax(1:(end-1)) - Ax(2:end);
V1y = Ay(1:(end-1)) - Ay(2:end);

V2x = Ax(1:(end-1)) - Bx_new(1:(end-1));
V2y = Ay(1:(end-1)) - By_new(1:(end-1));

V3x = Ax(1:(end-1)) - Bx_new(2:end);
V3y = Ay(1:(end-1)) - By_new(2:end);

Volumes_after = abs(V1x.*V3y - V3x.*V1y)/2 + abs(V2x.*V3y - V3x.*V2y)/2;
Lens_after = sqrt((Ax - Bx_new).^2 + (Ay - By_new).^2);

figure, hold on;
plot(Volumes_before, 'r');
plot(Volumes_after, 'b');
legend('Vol before', 'Vol after');

figure, hold on;
plot(Lens_before, 'r');
plot(Lens_after, 'b');
legend('Len before', 'Len after');

