% 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.    
%
%=========================================================================
% 
% Just type "DecomposeCurvature" without any arguments and things should work!
%
% This function performs a wavelet decomposition of a subject mesh. 
% It assumes the existence of WaveletBankBW500inv200.mat in the current
% directory. This .mat file is produced by first calling CreateLaplacianAnalysisSynthesisFilters.m 
% There are 5 input parameters you can play with.
%
% What to expect:
%           By default, there are 7 levels of filters produced by CreateLaplacianAnalysisSynthesisFilters.m 
%           The whole process should take only a few minutes.
%           At the end, you will see 7 figures, showing a multiresolution decomposition of the cortical curvature. 
%           You will also see 7 figures of the wavelet decomposition at each level plotted on the inflated cortical surface. 
%
% This function
%       1)  Calls the function PartiallyReconCurvature to do the wavelet decomposition 
%       2)  Save wavelet coefficients from each level onto "?h.curv.i", where i is the wavelet level.
%           Note that "?h.curv.i" is formed by passing the curvature through the i-th set of analysis filters. 
%       3)  Save each partially reconstructed multiresolution curvature into the "surf" folder of the subject
%           mesh as "?h.curv.recon.i" where i is the wavelet level.
%           Denote "?h.curv.syn.i" as the wavelet coefficients of level i "?h.wavelet.i" passing through the i-th synthesis filters. 
%           Now, ?h.curv.recon.j = \sum_{i=1}^j ?h.curv.syn.i. 
%           Thus "?h.curv.recon.j" is the partially reconstructed curvature up to level j
%       4)  Unlike the decomposition of the coordinate functions of the cortical surface in DecomposeSurface.m, 
%           the visual interpretation of the decomposition of curvature is much clearer. 
%           However, unlike the coordinate functions, we note that the mean curvature does not specify the entire surface. 
%       
% The general theories behind the overcomplete wavelet transform is based on the following papers (http://yeoyeo02.googlepages.com/publications).
% We note that the theories in the paper is far more general than the overcomplete wavelet codes we provide here.
% The papers discuss the construction of spherical steerable pyramid, which are overcomplete wavelet with oriented filters unlike the
% axisymmetric laplacian filters we use here. 
%
%           (a) On the Construction of Invertible Filter Banks on the 2-Sphere.  
%               B.T.T. Yeo, W. Ou and P. Golland. 
%               IEEE Transactions on Image Processing, 17(3):283--300, 2008
%
%           (b) Invertible Filter Banks on the 2-Sphere. 
%               B.T.T. Yeo, W. Ou and P. Golland. 
%               Proceedings of the IEEE International Conference on Image Processing (ICIP), 2161--2164, 2006
%
% The application of these overcomplete spherical wavelets to cortical surfaces is based on the following papers
% (http://yeoyeo02.googlepages.com/publications):
%
%           (c) Shape Analysis with Overcomplete Spherical Wavelets. 
%               B.T.T. Yeo, P. Yu, P.E. Grant, B. Fischl, P. Golland. 
%               Proceedings of the International Conference on Medical Image Computing and Computer Assisted Intervention (MICCAI), 
%               volume 5241 of LNCS, 468--476, 2008 
%
%           (d) Cortical Folding Development Study based on Over-complete Spherical Wavelets. 
%               P. Yu, B.T.T. Yeo, P.E. Grant, B. Fischl and P. Golland. 
%               Proceedings of the Workshop on Mathematical Methods in Biomedical Image Analysis, International Conference on Computer Vision, 2007
%

function DecomposeCurvature(hemi, SUBJECTS_DIR, subject, DISPLAY_RECONSTRUCTION, check_bool)

if(nargin < 5)
    subject_parms = CreateDefaultSubjectParms('lh');                     % left hemisphere ('lh') or right hemisphere ('rh')
    subject_parms.SUBJECTS_DIR = fullfile('..','..','example_surfaces'); % directory that contains subjects
    subject_parms.SUBJECT = 'OAS1_0001_MR1';                             % subject to be decomposed
    DISPLAY_RECONSTRUCTION = 1;                                          % DISPLAY_RECONSTRUCTION displays the partially reconstructed surfaces.
    check_bool = 0;                                                      % Makes checks and print stuff out (0 means no check).
else
    subject_parms = CreateDefaultSubjectParms(hemi);                           % left hemisphere ('lh') or right hemisphere ('rh')
    subject_parms.SUBJECTS_DIR = SUBJECTS_DIR;                           % directory that contains subjects
    subject_parms.SUBJECT = subject;                                     % subject in SUBJECTS_DIR
end

% Decomposes curvature.
subject_parms.data_filename_cell = {'curv'};
subject_parms.metric_surf_filename = 'inflated';
subject_parms.normalizeBool = 1; % 1 implies it will normalize the data. See MARS2_readSbjMesh.m for more details.

if(length(subject_parms.data_filename_cell) ~= 1)
   error('This function only decomposes one set of features'); 
end

% check to see if wavelets have already been created
if(ispc)
    check_dir = 'dir';
else
    check_dir = 'ls';
end
[a, b] = system([check_dir ' WaveletBankBW500inv200.mat']);

if(a ~= 0)
   error('WaveletBankBW500inv200.mat not found. Have you tried running CreateLaplacianAnalysisSynthesisFilters yet?');
else
    load('WaveletBankBW500inv200.mat');
end

% Perform Decomposition.
disp('Performing Wavelet Decomposition');
sbjMesh = feval(subject_parms.read_surface, subject_parms);
[recon_mat, analysis_mat] = PartiallyReconCurvature(sbjMesh, analysis_filters_my, synthesis_filters_my, check_bool);

% Saving the analysis coeffs.
disp('Saving Analysis Coefficients');
for j = 1:size(analysis_mat, 1)
    save_path = fullfile(subject_parms.SUBJECTS_DIR, subject_parms.SUBJECT, 'surf', [subject_parms.hemi '.' subject_parms.data_filename_cell{1} '.' num2str(j)]);
    write_curv(save_path, single(transpose(analysis_mat(j, :))), int32(size(sbjMesh.faces, 2))); 
end

% Saving the partially reconstructed curvature.
disp('Saving Partially Reconstructed Curvature');
for j = 1:size(recon_mat, 1)
    save_path = fullfile(subject_parms.SUBJECTS_DIR, subject_parms.SUBJECT, 'surf', [subject_parms.hemi '.' subject_parms.data_filename_cell{1} '.recon.' num2str(j)]);
    write_curv(save_path, single(transpose(recon_mat(j, :))), int32(size(sbjMesh.faces, 2))); 
end

% Display the partially reconstructed surfaces.
if(DISPLAY_RECONSTRUCTION)
    data_filename = subject_parms.data_filename_cell{1};
    for i = 1:size(recon_mat, 1)
        subject_parms.data_filename_cell = {[data_filename '.recon.' num2str(i)]};
        sbjMesh = feval(subject_parms.read_surface, subject_parms);
        sbjMesh.vertices = sbjMesh.metricVerts;
        s_handle = TrisurfMeshData(sbjMesh, sbjMesh.data, 1);
        shading interp;
        view(90, 0);
        camlight('headlight', 'infinite');
        view(270, 0);
        camlight('headlight', 'infinite');
        axis off;
        zoom(1.5);
        set(s_handle, 'DiffuseStrength', 1);
        set(s_handle, 'SpecularStrength', 0);
    end  
    
    for i = 1:size(analysis_mat, 1)
        subject_parms.data_filename_cell = {[data_filename '.' num2str(i)]};
        sbjMesh = feval(subject_parms.read_surface, subject_parms);
        sbjMesh.vertices = sbjMesh.metricVerts;

        s_handle = TrisurfMeshData(sbjMesh, sbjMesh.data, 1);
        shading interp;
        view(90, 0);
        camlight('headlight', 'infinite');
        view(270, 0);
        camlight('headlight', 'infinite');
        axis off;
        zoom(1.5);
        set(s_handle, 'DiffuseStrength', 1);
        set(s_handle, 'SpecularStrength', 0);
    end
end























