function [wfs beta]=WFS_COSINE(Y,x,k,sigma)
%
% function [wfs beta]=WFS_COSINE(Y,x,k,sigma)
%
% Construct the k-th degree weighted cosine  + sine series 
% representation on functions Y(x). This is equivalent 
% to performing heat kernel smoothing with bandwidth
% sigma. If sigma = 0, we have the usaul cosine Fourier series expansion.
%
% Y    : function value to smooth out
% x    : domain of function Y. dimension of Y match that of X.
% k    : degree of expansion
% sigma: amount of smoothing
%
% The technical detail and explanation on paramters are given in [1] and [2]
% 
% [1] Chung, M.K., Adluru, N., Lee, J.E., Lazar, M.,Lainhart, J.E., Alexander, A.L.
%     2010. Cosine series representation of 3D curves and its application to white 
%     matter fiber bundles in diffusion tensor imaging. Statistics and Its Interface. 
%     3:69-80. http://www.stat.wisc.edu/~mchung/papers/chung.2010.SII.pdf
%
% [2] Chung, M.K., Singh, V., Kim, P.T., Dalton, K.M., Davidson, R.J. 2009.
%     Topological characterization of signal in brain images using the min-max diagram.
%     12th International Conference on Medical Image Computing and Computer Assisted 
%     Intervention (MICCAI). Lecture Notes in Computer Science (LNCS). 5762:158-166.
%     http://www.stat.wisc.edu/~mchung/papers/miccai.2009.pdf 
%
% (C) Moo K. Chung 2008
%     mkchung@wisc.edu
%     University of Wisconsin-Madison
%
%     Histroy: 2008 created
%              2015 05/30 Input matrix data  

%n_vertex=length(Y);
n_vertex=size(Y,1);

% cosine basis
psi1=inline('sqrt(2)*cos(pi*l.*x)');
%psi2=inline('sqrt(2)*sin(pi*l.*x)');


%Y= psi*beta
% inv(psi)*Y = psi \ Y= beta

% design matrix
DX1=[];
for l=0:k
    psi=psi1(l,x);
    DX1=[DX1 psi];
end

% DX2=[];
% for l=1:k
%    psi=psi2(l,x);
%    DX2=[DX2 psi];
% end

%DX=[DX1 DX2];
DX=DX1;

beta=pinv(DX'*DX)*DX'*Y;

weight=diag(exp(-[0:k]'.^2*pi^2*sigma));
%weight = [exp(-[0:k]'.^2*pi^2*sigma) ;exp(-[1:k]'.^2*pi^2*sigma)];
wfs=DX*(weight*beta);


