%% Setup
close all
clear
clc

% Author: Christian Thode Larsen, christian.thode.larsen@gmail.com

% some environment dependent variables - configure as needed..
addpath '/opt/freesurfer/stable/matlab'

%% setup tests

basePath = '/home/cthla/Documents/phd/trunk/projects/bias/dev/data';

% 
% maskPath = '';
% maskS = MRIread(maskPath);

% ADNI Prepphase
dataPath = [basePath, '/ADNI_PrepPhase/007-S-0002-scan1.mgz'];
% NOTE! Prior is coregistered and resampled already
priorWMPath = [basePath, '/ADNI_PrepPhase/priors/007-S-0002-scan1-white.mgz'];
priorGMPath = [basePath, '/ADNI_PrepPhase/priors/007-S-0002-scan1-grey.mgz'];
priorCSFPath = [basePath, '/ADNI_PrepPhase/priors/007-S-0002-scan1-csf.mgz'];

dataS = MRIread(dataPath);
priorWM = MRIread(priorWMPath);
priorGM = MRIread(priorGMPath);
priorCSF = MRIread(priorCSFPath);

% prior is a 2D array, voxels x classes
prior = [priorWM.vol(:), priorGM.vol(:), priorCSF.vol(:)];

%runType = 'EM';
%runType = 'N3';
runType = 'sameVariance';
%runType = 'mixtureWithPrior';

showPlots = true;
%showPlots = false;

switch runType

    case 'EM'
        [correctedData, biasField, mask] = iic(dataS.vol, ...
            'stepSize', [dataS.xsize, dataS.ysize, dataS.zsize], ...
            'desiredStepSize', 4, ...
            'backgroundThresholdType', 'otzu', ...
            'numberOfGaussians', 6, ...
            'smoothingDistance', 50, ...
            'smoothingRegularization', 1e-7, ...            
            'showPlots', showPlots, ...
            'verbosity', 'debug');
    case 'N3'   
        [correctedData, biasField, mask] = iic(dataS.vol, ...
            'stepSize', [dataS.xsize, dataS.ysize, dataS.zsize], ...
            'desiredStepSize', 3, ...
            'backgroundThresholdType', 'otzu', ...
            'mimicN3', 'model', ...         
            'smoothingDistance', 50, ...
            'smoothingRegularization', 1e-7, ...               
            'showPlots', showPlots, ...
            'verbosity', 'debug');
    case 'sameVariance'
        [correctedData, biasField, mask] = iic(dataS.vol, ...
            'stepSize', [dataS.xsize, dataS.ysize, dataS.zsize], ...
            'desiredStepSize', 3, ...
            'backgroundThresholdType', 'otzu', ...
            'numberOfGaussians', 30, ...
            'sameVariance', true, ...
            'varianceScale', 1, ...
            'updateVariance', true, ...
            'equidistantMeans', true, ...
            'smoothingDistance', 50, ...
            'smoothingRegularization', 1e-7, ...              
            'showPlots', showPlots, ...
            'verbosity', 'debug');

    case 'mixtureWithPrior'
        [correctedData, biasField, mask] = iic(dataS.vol, ...
            'stepSize', [dataS.xsize, dataS.ysize, dataS.zsize], ...
            'desiredStepSize', 2, ...
            'prior', prior, ...
            'priorThreshhold', 0.99, ...
            'numberOfGaussians', [2,2,2,3], ...
            'smoothingDistance', 50, ...
            'smoothingRegularization', 1e-7, ...            
            'showPlots', showPlots, ...
            'verbosity', 'debug'); 
end
