function newY = MinArcLength(Mt,Mx,spikeI)
% minimize arc length of a cubic spline
% Mt: x-value of knots
% Mx: y-value of knots
% spikeI: index of knots to estimate new y-value
% ex: 
% [spmax, spmin, flag]=extrema(sig);
% Mx = spmax(:,2);
% Mt = spmax(:,1); 
% spikeI = 6;


% Code writer: Hui-Wen Yang 2019 @NTU and NCU
% Please cite:
% Hui-Wen Yang ; Shyh-Kang Jeng ; Hsu-Wen V. Young ; Chen Lin ; Yung-Hung Wang ; 
% Kun Hu ; Men-Tzung Lo. (2019). A Minimum Arclength Method for Removing Spikes 
% in Empirical Mode Decomposition. IEEE Access PP(99):1-1 · January 2019

% devide spline into 2 parts, 
% spline1 is (x0,y0), (x1,y1), ... (xi,0), (xi+1,yi+1), ... (xn,yn)
% spline2 is (x0, 0), (x1, 0), ... (xi,1), (xi+1,0), ... (xn,0)
% Target: find k s.t. spline1 + k*spline2 has the minimum arc length
x1 = Mx; x1(spikeI) = 0;
cs1 = spline(Mt,[0; x1; 0]);
x2 = zeros(size(Mx)); x2(spikeI) = 1;
cs2 = spline(Mt,[0; x2; 0]);

% % % to check the 2 spline equal
% dd=1:Mt(end);
% cs0 = spline(Mt,[0; Mx; 0]);
% sp0 = ppval(cs0,dd);
% sp1 = ppval(cs1,dd);
% sp2 = ppval(cs2,dd);
% sp12 = sp1 + Mx(spikeI)*sp2;
% figure; plot(sp0);
% hold all; plot(sp12)
% plot(Mt(spikeI),Mx(spikeI),'o')

a11 = cs1.coefs(spikeI-1,1);     b11 = cs1.coefs(spikeI-1,2);     c11 = cs1.coefs(spikeI-1,3);     d11 = cs1.coefs(spikeI-1,4);
a12 = cs2.coefs(spikeI-1,1);     b12 = cs2.coefs(spikeI-1,2);     c12 = cs2.coefs(spikeI-1,3);     d12 = cs2.coefs(spikeI-1,4);
a21 = cs1.coefs(spikeI,1);     b21 = cs1.coefs(spikeI,2);     c21 = cs1.coefs(spikeI,3);     d21 = cs1.coefs(spikeI,4);
a22 = cs2.coefs(spikeI,1);     b22 = cs2.coefs(spikeI,2);     c22 = cs2.coefs(spikeI,3);     d22 = cs2.coefs(spikeI,4);
% Arc Length = intergral { sqrt(1+dy/dx) } dx
syms x k
pp1 = sqrt((3*(a11+k*a12)*x^2 + 2*(b11+k*b12)*x + c11+k*c12)^2+1);
pp2 = sqrt((3*(a21+k*a22)*x^2 + 2*(b21+k*b22)*x + c21+k*c22)^2+1);

ArcLength = int(pp1,x,Mt(spikeI-1),Mt(spikeI)) + int(pp2,x,Mt(spikeI),Mt(spikeI+1));
g = diff(ArcLength,k);
S = solve(g);
% Sm = subs(ArcLength, k, S);
% Mx(6) = S;
newY = S;
% figure; plot(tsig,sig);
% hold on;
% plot(tsig(Mt),Mx,'o')
% plot(tsig(spmin(:,1)),spmin(:,2),'o')
%% 


% f = @(x) sqrt(1+dydx(x));
% arclength = intergral(f, t_a, t_b);
