% Contact ythomas@csail.mit.edu or msabuncu@csail.mit.edu for bugs or questions 
%
%=========================================================================
%
%  Copyright (c) 2008 Thomas Yeo and Mert Sabuncu
%  All rights reserved.
%
%Redistribution and use in source and binary forms, with or without
%modification, are permitted provided that the following conditions are met:
%
%    * Redistributions of source code must retain the above copyright notice,
%      this list of conditions and the following disclaimer.
%
%    * Redistributions in binary form must reproduce the above copyright notice,
%      this list of conditions and the following disclaimer in the documentation
%      and/or other materials provided with the distribution.
%
%    * Neither the names of the copyright holders nor the names of future
%      contributors may be used to endorse or promote products derived from this
%      software without specific prior written permission.
%
%THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
%ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
%WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
%DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
%ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
%(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
%LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
%ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
%(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
%SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.    
%
%=========================================================================
%
% This function decomposes the feature a closed 2D subject mesh in a multiresolution fashion (MARS_sbjMesh.data).
% The bulk of the work is done by WaveletTransformCurvature.m
%
% WaveletTransformCurvature returns the outputs of the analysis filters (wc) 
% and the outputs of the synthesis filters (rc)
%
% Since wc and rc are expressed in spherical harmonic coefficients, the function goes on to
% perform an inverse spherical harmonic transform and returns a series of partially reconstructed features 
% and the wavelet coefficients of the features at each spatial level.

function [recon_mat, analysis_mat] = PartiallyReconCurvature(MARS_sbjMesh, analysis_filters, synthesis_filters, check_bool)

% [recon_cell] = PartiallyReconCurvature(MARS_sbjMesh, analysis_filters, synthesis_filters, check_bool)

if(size(analysis_filters, 2) ~= size(synthesis_filters, 2))
   error('We expect analysis and synthesis filters to be of same bandwidth');
end

if(size(analysis_filters, 1) ~= size(synthesis_filters, 1))
   error('We expect the number of analysis and synthesis filters to be the same');
end

% wavelet transformation
num_filters = size(analysis_filters, 1);
[wc, rc] = WaveletTransformCurvature(MARS_sbjMesh, analysis_filters, synthesis_filters, check_bool);

% reconstruct brain
coeffs_mat = zeros(size(rc));

% Compute theta,phi of lat-lon grid
BW = 1/2 *(-1 + sqrt(1+8*size(analysis_filters, 2)));
[phi, theta] = sphgrid(BW*2);
phi = [phi 2*pi*ones(size(phi, 1), 1)];
phi = [phi(1,:)   ; phi;  phi(1,:)]; 
theta = [theta theta(:, 1)];
theta = [zeros(1, size(theta, 2)) ; theta; pi*ones(1, size(theta, 2))];

% Compute theta,phi of cortical surface mesh (which had been spherically parameterized)
radius = sqrt(sum(MARS_sbjMesh.vertices.^2, 1));
offset_vec = zeros(1, size(MARS_sbjMesh.vertices, 2)); 
offset_vec(abs(MARS_sbjMesh.vertices(3, :) - 100) < 1e-5) = 1e-5; % These offsets arise due to numerical issues so that we won't get cos or sin larger than 1
offset_vec(abs(MARS_sbjMesh.vertices(3, :) + 100) < 1e-5) = -1e-5;% Or else I will get complex numbers for theta and phi.
theta_vec = mod(single(acos(MARS_sbjMesh.vertices(3,:)./radius)) + offset_vec, single(pi));
phi_vec = mod(single(atan2(MARS_sbjMesh.vertices(2,:),MARS_sbjMesh.vertices(1,:))), single(2*pi));


analysis_mat = zeros(num_filters, size(MARS_sbjMesh.vertices, 2));
for i = 1:num_filters
    
    disp(['Performing Inverse Spherical Harmonic Transform of Wavelet Coefficients level ' num2str(i) ' followed by interpolation onto original mesh points']);
    
    mat = ConvertMyFormat2S2kitFormat(transpose(wc(i, :)));
    mesh_s2kit_coeff = s2ifst(mat);
    mesh_s2kit_coeff = [mesh_s2kit_coeff mesh_s2kit_coeff(:, 1)];
    mesh_s2kit_coeff = [repmat(mean(mesh_s2kit_coeff(1,:)), 1, size(mesh_s2kit_coeff, 2)); mesh_s2kit_coeff; repmat(mean(mesh_s2kit_coeff(end,:)), 1, size(mesh_s2kit_coeff, 2))];
    analysis_mat(i, :) = MARS_interp2(phi, theta, real(mesh_s2kit_coeff), phi_vec, theta_vec, '*linear');    
end

recon_mat = zeros(num_filters, size(MARS_sbjMesh.vertices, 2));
for i = 1:num_filters
    
    disp(['Performing Inverse Spherical Harmonic Transform of Reconstruction level ' num2str(i) ' followed by interpolation onto original mesh points']);
    coeffs_mat(i, :) = sum(rc(1:i, :), 1);
    
    mat = ConvertMyFormat2S2kitFormat(transpose(coeffs_mat(i, :)));
    mesh_s2kit_recon = s2ifst(mat);
    mesh_s2kit_recon = [mesh_s2kit_recon mesh_s2kit_recon(:, 1)];
    mesh_s2kit_recon = [repmat(mean(mesh_s2kit_recon(1,:)), 1, size(mesh_s2kit_recon, 2)); mesh_s2kit_recon; repmat(mean(mesh_s2kit_recon(end,:)), 1, size(mesh_s2kit_recon, 2))];
    recon_mat(i, :) = MARS_interp2(phi, theta, real(mesh_s2kit_recon), phi_vec, theta_vec, '*linear');
end






