function [averagecompo,averagecompovf]=dibujopca(study_fourier,sphere, k,sigma, factorm,desvcoi,coi,r);
% 
% study_fourier  : A structured array of Fourier coefficients. 
%                        Example: The Fourier coefficinets of the 10th-subject is given by
%
%                        > study_fourier(10,:)
%                             ans =
%                                   x: [43x85 double]
%                                   y: [43x85 double]
%                                   z: [43x85 double]
%                
% sphere:  A unit spherical mesh where the coordinates will be sampled.
% sigma: bandwidth of smoothing
% k : degree. For amygdala, k=15, cortical surface k = 48 or 52. 
%
% factorm: es el vector de factores por el que multiplicaré la desviación típica de cada
% compoenente, para su representacion: factor multiplying sd each component
%
% desvcoi: desviación típica de la componente que represento: sd of the component
%
% coi, es el vector de coeficientes de la componente: coefficient vector of the component
%
% r: donde colocar los coeficientes en base de Fourier: where the coef. are in the base
%
% Devuelve, los coeficientes y hace la representacion para ese factor: devuelve vértices y
% caras: it returns the coefficients and the surface for that component and factor


n_surfaces=size(study_fourier,1);

% chop off up to k-th degree
for i=1:n_surfaces
    temp=study_fourier(i);
    study_fourier(i).x=temp.x(1:k+1,1:k*2+1);
    study_fourier(i).y=temp.y(1:k+1,1:k*2+1);
    study_fourier(i).z=temp.z(1:k+1,1:k*2+1);
end;


left_x=0;left_y=0;left_z=0;

for i=1:n_surfaces
    left=study_fourier(i);
    left_x = left_x + left.x;
    left_y = left_y + left.y;
    left_z = left_z + left.z;
end;

left_x=left_x/n_surfaces; 
left_y=left_y/n_surfaces; 
left_z=left_z/n_surfaces;

fourier_coeff.x=left_x;
fourier_coeff.y=left_y;
fourier_coeff.z=left_z;

cx=convenmatrix(coi(1:256),r);

cy=convenmatrix(coi(257:512),r);
cz=convenmatrix(coi(513:768),r);


averagecompo.x= fourier_coeff.x + factorm(1)*desvcoi*cx;
averagecompo.y= fourier_coeff.y + factorm(2)*desvcoi*cy;
averagecompo.z= fourier_coeff.z + factorm(3)*desvcoi*cz;

%averagecompovf=[];


averagecompovf=SPHARMrepresent2(sphere,averagecompo,k,sigma);





%----------------------------------------------------------------------
function cmatriz=convenmatrix(vectc,r)

cmatriz=zeros(16,31);
cmatriz=cmatriz(:);
cmatriz(r)=vectc;

cmatriz=reshape(cmatriz,16,31);


%--------------------------------------------------------------------------
function SPHARMsurf=SPHARMrepresent2(surf,fourier_coeff,k,sigma)

n_vertex = size(surf.vertices,1);

% determining the number of subjects
if ndims(fourier_coeff.x)==3
    n_subject=size(fourier_coeff.x,3);
else
    n_subject=1;
end;

[theta varphi]=EULERangles(surf);

% Initialization for coordinates
xcoord=zeros(n_vertex,1);
ycoord=zeros(n_vertex,1);
zcoord=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=fourier_coeff.x(1,1);
xcoord=Y*betal;

betal=fourier_coeff.y(1,1);
ycoord=Y*betal;

betal=fourier_coeff.z(1,1);
zcoord=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=fourier_coeff.x(l+1,1:2*l+1);
    xsmooth=Y*betal';
    xcoord=xcoord+ exp(-l*(l+1)*sigma)*xsmooth;

    betal=fourier_coeff.y(l+1,1:2*l+1);
    ysmooth=Y*betal';
    ycoord=ycoord+ exp(-l*(l+1)*sigma)*ysmooth;

    betal=fourier_coeff.z(l+1,1:2*l+1);
    zsmooth=Y*betal';
    zcoord=zcoord+ exp(-l*(l+1)*sigma)*zsmooth;

end;

%output the results in a proper shape

temp=[xcoord ycoord zcoord];
SPHARMsurf.vertices=squeeze(reshape(temp,n_vertex,3));
SPHARMsurf.faces=surf.faces;

%---------------------------------------------------------------
function [theta,varphi]=EULERangles(surf);

n_vertex=size(surf.vertices,1);
c=mean(surf.vertices);  %mass center
surf.vertices=surf.vertices-kron(ones(n_vertex,1),c);  % translation

[theta,varphi,r] = cart2sph(surf.vertices(:,1),surf.vertices(:,2),surf.vertices(:,3));

% MATLAB coordinate systems are different from the convention used in the
% TMI paper.
temp = theta;
theta = pi/2 - varphi;
varphi = pi + temp;

%figure_wire(surf,'yellow')

%-----------------------------------------------------------------
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;




