function [square, xcoord]=SPHARM2square(coeff,k,sigma)
%
%function square=SPHARM2square(coeff,k,sigma)
%
% Project the weighted spherical harmonic representation onto a square
%
%
%
% coeff   : The estimated SPHARM-coefficients 
%     k   : degree of expansion
% sigma   : bandwidth
%
%
% (C) Moo K. Chung 2008-2010
%     Department of Biostatistics and Medical Informatics
%     Waisman Laboratory for Brain Imaging
%     University of Wisconsin-Maison
%  
% email://mkchung@wisc.edu
% http://www.stat.wisc.edu/softwares/weighted-SPHARM/weighted-SPHARM.html
%
% If you use this code, please reference the following paper. 
% You need to read the paper to understand the notations and the algorithm.
%
% Chung, M.K., Dalton, K.M., Shen, L., L., Evans, A.C., Davidson, R.J. 2007. 
% Weighted Fourier series representation and its application to quantifying 
% the amount of gray matter. IEEE Transactions on Medical Imaging, 26:566-581.
%
% July 6, 2007 created 
% Jan  9, 2008 modified; July 19, 2010: simplified
%----------------------------------------------------------------------------------------

theta_d=0:0.01:pi;
varphi_d=0:0.01:2*pi;
m=length(theta_d);
n=length(varphi_d);
n_vertex=m*n;

theta=kron(ones(1,n),theta_d');
theta=reshape(theta,m*n,1);

varphi=kron(ones(m,1),varphi_d);
varphi=reshape(varphi,m*n,1);

% Initialization for coordinates

n_vertex=length(theta);
xcoord=zeros(n_vertex,1);

% 0-th degree construction

% real(Y_l) gives negative harmonics imag(Y_l) gives positive harmonics
% reading Y_lm and constructing design matrix. See the paper.
Y=Y_l(0,theta',varphi')';

betal=coeff(1,1);
xcoord=Y*betal;


% Iteratively add each degree up to the k-th degree.    

for l=1:k

    Y=Y_l(l,theta',varphi')';
    % real(Y_l) gives negative harmonics imag(Y_l) gives positive harmonics
    Y=[real(Y)   imag(Y(:,2:(l+1)))];
    
    betal=coeff(l+1,1:2*l+1);
    xsmooth=Y*betal'; 
    xcoord=xcoord+ exp(-l*(l+1)*sigma)*xsmooth; 
    
end;

%output the results in a proper shape

%SPHARMsurf.vertices=squeeze(reshape(temp,n_vertex,3));

square=xcoord;

square=reshape(square,m,n);


%-----------------------------------------------------------------
function Y_l=Y_l(l,theta,varphi)
% computes spherical harmonics of degree l.
sz=length(theta); 

m=0:l;
CLM=[];
exp_i_m=[];
sign_m=[];
SIGNM=[];
Pn=[];

for k = 0:(2*l)
    fact(k+1) = factorial(k);
end
clm = sqrt(((2*l+1)/(2*pi))*(fact(l-abs(m)+1)./fact(l+abs(m)+1)));
CLM=kron(ones(1,sz),clm');

for k = 0:l
    exp_i_m(k+1,:)= exp(i*k*varphi);
    sign_m(k+1) = (-1)^k;
end
exp_i_m(1,:)=exp_i_m(1,:)/sqrt(2);

SIGNM=kron(ones(1,sz),sign_m');
Pn=legendre(l,cos(theta));
Y_l=CLM.*SIGNM.*Pn.*exp_i_m;