function varargout = norm_fluid_auto(varargin)
% CG_VBM_OPTIMIZED
% optimized processing of voxel based morphometry
%
% It is recommended to use your own gray matter template:
% First, run the script cg_create_template with all raw images from your sample.
%
% inputs
% V                 : volumes to normalize       (spm volume structure)
% VGO               : template for normalisation (spm volume structure)
% graphicmode       : 0|1 (no graphics | graphics)
% estimate          : 0|1 (average previous estimates | make new estimates & average)
% modulation        : 0|1 (no|yes)
% optimal           : 0|1 (standard VBM | optimized VBM)
% NRMsourceweight   : 0|1 (no|yes)
% NRMtemplateweight : 0|1 (no|yes)
% SEGsourceweight   : 0|1 (no|yes)
% SEGtemplateweight : 0|1 (no|yes)
% ApplyMask         : 0|1 (no|yes)
% write_seg         : 0|1 (no|yes) write_seg=1 saves segmented GM, WM, and CSF images to disk
%
%
% output files:
%---------------
%	w*.img		warped raw image
%	Gw*.img:	gray matter of warped raw image
%	Ww*.img:	white matter of warped raw image
%	Cw*.img:	CSF of warped raw image
%	m[GWC]w*.img:	modulated segmented image of warped raw image
%	avgT1.img	average of all warped T1 images
% 	vbm_volumes.txt	absolute volumes of gray/white matter and csf
%
% modulation:
%------------
% This step is intended to preserve the total amount of tissue in the images.
% Areas that are expanded during warping are correspondingly reduced
% in intensity. John Ashburner recommends this step if you are more
% interested in volume changes than differences in concentration
% (or density). Modulates images will be additionally saved indicated
% by a leading 'm'. Nevertheless this postprocessing is optional and you can
% try both, with and without modulation. If you choose this post-
% processign step, modulated and non-modulated results will be saved.
%
% use of cutomized prior images
% -----------------------------
% Depending on your dataset it is also useful to use cutomized prior
% images resulting from cg_create_template.
% The use of own priors is recommended if the gray/white matter
% distribution of your sample differs from a control population
% (e.g. children, Alzheimer disease...).
%
% calculation of absolute gray/white matter and csf volumes
% ---------------------------------------------------------
% A file "vbm_volumes.txt" is saved in the actual directory with the
% absolute volumes of gray/white matter and CSF (in ml). If this file
% exist the new values will be appended. This file can be opened with
% Excel (or gnumeric) as text-file with tabs as delimiters.
%_______________________________________________________________________
% @(#)cg_vbm_optimized.m	1.18 Christian Gaser 2004/03/04
% modified script of John Ashburner
%  spm_defaults
[V,matfiles,W,VG0,priors,graphicmode,modulation,optimal,NRMsourceweight,...
    NRMtemplateweight,SEGsourceweight, SEGtemplateweight, ApplyMask, ...
    max_seg, writeseg,writeanat,writeepi,bmask,voxsize] = deal(varargin{1:nargin});
disp('normalizing ... (norm_fluid_auto.m)')


writeaverage=1;
own_prior=1;
global defaults
dseg				= defaults.segment;
dseg.write.wrt_cor		= 0;		% don't write bias corrected image
%dseg.estimate.affreg.weight	= fullfile(spm('Dir'),'apriori','brainmask.mnc');	% always use brainmask
dseg.write.cleanup   		= 1;		% clean up partitions
dseg.estimate.affreg.sourceweight = '';
VWF='';
VWG='';

dnrm				= defaults.normalise;
if ~isempty(voxsize)
    switch length(voxsize)
        case 1
            dnrm.write.vox			= [voxsize voxsize voxsize];%[.859 .859 .859];%[1 1 1];
        case 3
            dnrm.write.vox			= voxsize;%[.859 .859 .859];%[1 1 1];
    otherwise
    error('voxsize must be length 1 or 3!!!')
    end
end
%dnrm.write.bb(1,3)		= -70;		% enlarge z-axis to include more cerebellum
dnrm.write.bb = [Inf Inf Inf; Inf Inf Inf];
%dnrm.write.vox			= [1 1 1];
dnrm.estimate.graphics		= 1;		% Crashes otherwise
if ~exist('V','var')
    V   = spm_vol(spm_get(Inf,'*.IMAGE','Select raw (not normalized) images'));
end
if ~exist('matfiles','var')
    matfiles   = {};%spm_get([0 length(V)],'*_sn.mat',{'Select *_sn.mat files'});
end
if ~exist('W','var')
    for i=1:length(V)
        [pth,nam,ext] = fileparts(V(i).fname);
        str=['Select EPI images for ',nam];
        fltr=['*',nam(1:4),'*.IMAGE'];
        WW  = spm_get(Inf,{fltr},{str});
        if ~isempty(WW)
            W(i,:)   = spm_vol(WW);
        end
    end
end
if ~exist('VG0','var')
    VG0 = spm_vol(spm_get(1,'*T1*.IMAGE','Select T1 template for first normalization'));
end
if ~exist('max_seg','var')

    max_seg    = spm_input('Calculate',1,'m','gray matter|gray+white matter|all|all + intermediate images',[1 2 3 4],1);
end
if ~exist('modulation','var')

    modulation = spm_input('Additional modulation?','+1','yes|no',[1 0],2);
end
if ~exist('own_prior','var')
    own_prior  = spm_input('Use own prior images?','+1','yes|no',[1 0],2);
end
if ~exist('priors','var')

    if own_prior
        dseg.estimate.priors = spm_get(3,'*.IMAGE','Select gray/white and CSF prior images');
    end
else
    dseg.estimate.priors =priors;
end

% change default parameters for (non)linear normalization
%dnrm.estimate.cutoff = get_cutoff(dnrm.estimate.cutoff);
%if dnrm.estimate.cutoff ~= Inf
%    dnrm.estimate.reg    = get_reg(dnrm.estimate.reg);
%    dnrm.estimate.nits   = get_nits(dnrm.estimate.nits);
%end
if ~exist('NRMsourceweight','var')
    NRMsourceweight   = get_weight(0,'NRM source');
end
if ~exist('NRMtemplateweight','var')
    NRMtemplateweight = get_weight(0,'NRM template');
end
if ~exist('SEGsourceweight','var')
    SEGsourceweight   = get_weight(0,'SEG source');
end
if ~exist('SEGtemplateweight','var')
    SEGtemplateweight = get_weight(0,'SEG template');
end
if ~exist('ApplyMask','var')
    ApplyMask = get_ApplyMask(0);
end
if NRMtemplateweight
    VWG = fullfile(spm('Dir'),'apriori','brainmask.mnc');%
else
    VWG = '';
end
if SEGtemplateweight
    dseg.estimate.affreg.weight= fullfile(spm('Dir'),'apriori','brainmask.mnc');%
else
    dseg.estimate.affreg.weight= '';
end
%if ~exist('VG1','var')
%    VG1 = spm_vol(spm_get(1,'*gray*.IMAGE','Select own prior (gray) image for normalization'));
%end
if optimal==0
    VG1 = VG0;
elseif optimal==1
    VG1 = spm_vol(dseg.estimate.priors(1,:));
elseif optimal==2
    VG1 = spm_vol(dseg.estimate.priors(2,:))
end

% create graphics window
fg=spm_figure('FindWin','Graphics');
if isempty(fg), fg=spm_figure('Create','Graphics'); end

if exist('reuse_anat','var') && reuse_anat
else
    % save gray/white matter and csf volumes in log-file
    logname = 'vbm_volumes.txt';
    if exist(logname,'file');
        fid = fopen(logname,'a');
    else
        fid = fopen(logname,'w');
        fprintf(fid,'Volume [ml]\tgray\twhite\tcsf');
    end
    save NFAparams
    disp('saving parameters in NFAparams.mat');

    n = (0);
    l=length(V);
    for i=1:l
        disp(['Working on image ',num2str(i),'/',num2str(l)])
        workdir=pwd;
        save norm_fluid_auto_restart_params;
        %%%
        if (NRMsourceweight | SEGsourceweight | ApplyMask)
            [pth,nam,ext,patid] = GetPatientID(V(i).fname);
        else
            [pth,nam,ext] = fileparts(V(i).fname);
        end
        matname0 = fullfile(pth,[nam '_sn.mat']);
        if exist(matname0,'file')
            m0=load(matname0);
        end

        % get the patient ID from the filename
        % apply manually edited brain mask to improve the normalisation.
        if NRMsourceweight
            VWF = fullfile(pth,'edits','bin', ['MSK_',patid,'.img']);
            if ApplyMask, msk=spm_vol(VWF); end
        end

        % in mls_normalise_DEV:
        % VWG       - template weighting image
        % VWF       - source weighting image
        %VWF = dnrm.estimate.weight;
        if  SEGsourceweight
            dseg.estimate.affreg.sourceweight =fullfile(pth,'edits','bin', ['MSK_',patid,'.img']);
            if ApplyMask, msk = spm_vol(dseg.estimate.affreg.sourceweight); end
        end
        if exist('matfiles','var') & ~isempty(matfiles)
            matname1=matname0;
            try
                matname0 = matfiles{i,:};%pth,[nam '_sn.mat']);
            catch
                matname0 = matfiles{i};
            end

            m0=load(matname0);
            matname1 = fullfile(pth,[nam '_sn2.mat']);
        end
        if optimal
            % Estimate spatial normalisation parameters
            % by matching GM with a GM template
            disp(['segmenting ',nam,ext])
            if exist('m0','var') & isfield(m0,'M')
                disp(['using initial M estimate from ',matname0])
                sflags=dseg;
                rerun=1;
                while rerun<4 & rerun
                    %disp(['rerun3 = ',num2str(rerun)])
                    [VT,V(i),rerun,sflags] = spm_segment_mls2(V(i),m0.M,sflags,rerun);
                end
                if rerun==4
                    rerun=0;
                    [VT,V(i),rerun,sflags] = spm_segment_mls2(V(i),m0.M,sflags,rerun);
                end
            else
                sflags=dseg;
                rerun=1;
                while rerun<3 & rerun
                    %disp(['rerun4 = ',num2str(rerun)])
                    [VT,V(i),rerun,sflags] = spm_segment_mls2(V(i),VG0,sflags,rerun);
                end
                if rerun>=3
                    rerun=0;
                    [VT,V(i),rerun,sflags] = spm_segment_mls2(V(i),VG0,sflags,rerun);
                end
            end
            if optimal==1 || max_seg==4
                VT(1).fname = fullfile(pth,['tmpGM_'  nam ext]);
                try
                    spm_write_vol(VT(1),double(VT(1).dat/255));
                catch
                    rerun = 0;
                    save
                    [VT,V(i),rerun,sflags] = spm_segment_mls2(V(i),VG0,sflags,rerun);
                    save
                    %VT(1).fname = fullfile(pth,['tmpGM_'  nam ext]);
                    spm_write_vol(VT(1),double(VT(1).dat/255));
                    save
                end
            end
            if optimal==2 || max_seg==4
                VT(2).fname = fullfile(pth,['tmpWM_'  nam ext]);
                spm_write_vol(VT(2),double(VT(2).dat/255));
            end
            if optimal==1,
                VTname=VT(1).fname;             
            elseif optimal==2,  
                VTname=VT(2).fname; 
            end
        else
            VTname=V(i).fname;
            VG1 = VG0;
        end
        clear VT
        VT=spm_vol(VTname);
        % if (ApplyMask)
        %     % apply manually edited brain mask to improve the segmentation
        %     VTo  = ApplyBrainMask(VT,msk);
        %     clear VT
        %     VT = VTo;
        %     clear VTo
        % end
        disp(['normalizing ',VTname])
        if  exist('matname0','var') & exist(matname0,'file')
            prm = mls_normalise_DEV_2(VG1,VT,matname1,VWG,VWF,dnrm.estimate,matname0);
        elseif  exist('matname0','var') & ~exist(matname0,'file')
            prm = mls_normalise_DEV_2(VG1,VT,matname0,VWG,VWF,dnrm.estimate);
        else
            error('matname0 not defined')
        end
        clear VT
        if optimal,
            % remove temporary GM / WM image
            %%mls_cleanup({VTname});
        end
        % Create a spatially normalised version of the original image & brainmask
        disp(['writing normalized Anatomical image ...']);
        if isempty(voxsize)
            dnrm.write.vox = abs(GetVoxSize(V(i)));
        end
        VN  = spm_write_sn(V(i),prm,dnrm.write);
        VN.fname = fullfile(pth,['w_' nam ext]);
        if writeanat
            spm_write_vol(VN,VN.dat);
        end
        if writeaverage
            try
                % make average of warped images
                n = n + double(VN.dat);
            catch
                writeaverage=0;
                warning('couldn''t average images')
            end
        end
        %    if modulation & writeanat
        %        disp(['writing modulated normalized Anatomical image ...']);
        %        spm_write_sn(VN,prm,'modulate');
        %    end
        if writeepi
            if exist('W','var')
                disp(['writing normalized EPI / PET / other images ...']);
                if iscell(W)
                    k=size(W,2);
                    for j=1:k
                        %	  disp(['Saving Image ', num2str(j), '/',num2str(k)]);
                        % Print progress information in command window
                        str   = sprintf('Saving Image %3d/%-3d ',j,k);
                        fprintf('\r%-40s: %20s',str,' ')
                        try www=W{i,j};
                        catch www=W{i};
                        end
                        [Wpth,Wnam,Wext] = fileparts(www.fname);
                        if isempty(voxsize)
                            dnrm.write.vox = abs(GetVoxSize(www));
                        end
                        WN = spm_write_sn(www,prm,dnrm.write);
                        WN.fname = fullfile(Wpth,['w_' Wnam Wext]);
                        spm_write_vol(WN,WN.dat);
                        %%delname{j} =  fullfile(Wpth,[Wnam Wext]);
                        clear WN
                    end
                elseif isstruct(W)
                    [k1,k2]=size(W);
                    if k1==length(V)
                        k=k2;
                        for j=1:k
                            %    disp(['Saving Image ', num2str(j), '/',num2str(k)]);
                            % Print progress information in command window
                            str   = sprintf('Saving Image %3d/%-3d ',j,k);
                            fprintf('\r%-40s: %20s',str,' ')

                            [Wpth,Wnam,Wext] = fileparts(W(i,j).fname);
                            if isempty(voxsize)
                                dnrm.write.vox = abs(GetVoxSize(W(i,j)));
                            end
                            WN = spm_write_sn(W(i,j),prm,dnrm.write);
                            WN.fname = fullfile(Wpth,['w_' Wnam Wext]);
                            %%delname{j} =  fullfile(Wpth,[Wnam Wext]);
                            spm_write_vol(WN,WN.dat);
                            clear WN
                        end
                    elseif k2==length(V)
                        k=k1;
                        for j=1:k
                            %     disp(['Saving Image ', num2str(j), '/',num2str(k)]);
                            % Print progress information in command window
                            str   = sprintf('Saving Image %3d/%-3d ',j,k);
                            fprintf('\r%-40s: %20s',str,' ')

                            [Wpth,Wnam,Wext] = fileparts(W(j,i).fname);
                            if isempty(voxsize)
                                dnrm.write.vox = abs(GetVoxSize(W(j,i)));
                            end
                            WN = spm_write_sn(W(j,i),prm,dnrm.write);
                            WN.fname = fullfile(Wpth,['w_' Wnam Wext]);
                            spm_write_vol(WN,WN.dat);
                            %%delname{j} =  fullfile(Wpth,[Wnam Wext]);
                            clear WN
                        end
                    end
%                    %%%%mls_cleanup(%%delname);
                    %%clear %%delname
                end
            end
            fprintf('%s%20s',repmat(sprintf('\b'),1,20),'... done',sprintf('\n'))%-#
        end
        if ApplyMask
            % save the mask image
            disp(['writing normalized Mask image ...']);
            BrainMaskN  = spm_write_sn(msk,prm,dnrm.write);
            BrainMaskN.fname = fullfile(pth,['w_BrainMask' nam ext]);
            foo = BrainMaskN.fname;
            spm_write_vol(BrainMaskN, BrainMaskN.dat);
            clear BrainMaskN;
            BrainMaskN = spm_vol(foo);
        end
        % Segment the spatially normalised image
        disp(['segmenting ',VN.fname])
        sflags=dseg;
        rerun=1;
        while rerun<4 & rerun
            %disp(['rerun2 = ',num2str(rerun)])
            [VT,VN,rerun,sflags]  = spm_segment_mls2(VN,eye(4),sflags,rerun);
        end
        VT(1).fname = fullfile(pth,['tmpGM2_'  nam ext]);
        VT(2).fname = fullfile(pth,['tmpWM2_'  nam ext]);
        spm_write_vol(VT(1),double(VT(1).dat/255));
        spm_write_vol(VT(2),double(VT(2).dat/255));
        VT1name=VT(1).fname;
        VT2name=VT(2).fname;
        VT1=spm_vol(VT1name);
        VT2=spm_vol(VT2name);
        VT_3=VT(3);
        if (ApplyMask)
            disp('applying normalized brain mask to improve segmentation')
            VT_1  = ApplyBrainMask(VT1,BrainMaskN);
            VT_2  = ApplyBrainMask(VT2,BrainMaskN);
            mfname{1} = fullfile(pth,['btmpGM2_'  nam ext]);
            mfname{2} = fullfile(pth,['btmpWM2_'  nam ext]);
            VT_1.dat = uint8(get_voxel_values3d(VT_1,VT_1.dim(3),0));
            VT_2.dat = uint8(get_voxel_values3d(VT_2,VT_2.dim(3),0));
            % VT_3.dat = get_voxel_values3d(VT_3,VT_3.dim(3),0);
        else
            VT_1 = VT(1);
            VT_2 = VT(2);
        end
        VT_1.fname = fullfile(pth,['Gw_'  nam ext]);
        VT_2.fname = fullfile(pth,['Ww_'  nam ext]);
        VT_3.fname = fullfile(pth,['Cw_'  nam ext]);

        % CSF is better estimated using gray and white matter values
        t1 = ((double(VT_1.dat) + double(VT_2.dat) + double(VT_3.dat)) > 0.5*255);
        t2 = (1 - (double(VT_1.dat) + double(VT_2.dat))/255);
        VT_3.dat = uint8(round((t1.*t2)*255));

        % calculate volumes of GM, WM and CSF
        factor = abs(det(VT_1.mat(1:3,1:3)))/1000;
        gray  = factor*sum(VT_1.dat(:))/255;
        white = factor*sum(VT_2.dat(:))/255;
        csf   = factor*sum(VT_3.dat(:))/255;

        % save volumes in log-file
        fprintf(fid,'\n%s\t%5.3f\t%5.3f\t%5.3f\t%5.3f',nam,gray,white,csf);

        %  VT_1.dat = uint8(round(double(VT_1.dat)));
        %  VT_2.dat = uint8(round(double(VT_2.dat)));

        if size(VT_1.pinfo,1)>2, VT_1.pinfo(3:end)=[]; end
        if size(VT_2.pinfo,1)>2, VT_2.pinfo(3:end)=[]; end
        if size(VT_3.pinfo,1)>2, VT_3.pinfo(3:end)=[]; end
        if writeseg
            if modulation,
                disp(['writing modulated segmented images ...']);
                if max_seg>0, spm_write_sn(VT_1,prm,'modulate');end
                if max_seg>1, spm_write_sn(VT_2,prm,'modulate');end
                if max_seg>2, spm_write_sn(VT_3,prm,'modulate');end
            end
            disp(['writing unmodulated segmented images ...']);
            if max_seg>0,spm_write_vol(VT_1,double(VT_1.dat)/255);end
            if max_seg>1,spm_write_vol(VT_2,double(VT_2.dat)/255);end
            if max_seg>2,spm_write_vol(VT_3,double(VT_3.dat)/255);end
        end
        if max_seg<=3
            %%%%mls_cleanup({VT1name,VT2name});
        else
            disp(['saved the following images for pet PVE:'])
            disp(VT1name)
            disp(VT2name)
            PETpve = {VT1name,VT2name}
        end
        if exist('mfname','var')
            %%%%mls_cleanup(mfname);
        end
        clear VT
    end
end
%
%%%

fclose(fid);
if writeaverage
    % save average image
    n = n/length(V);
    VN.fname   = 'avgT1.img';
    VN.descrip = 'Average of warped T1 images';
    spm_write_vol(VN,n);
else
    warning('could not write average T1 image... ')
    disp('voxel sizes are probably not all the same')
end
try
    delete norm_fluid_auto_restart_params.mat;
end
if max_seg==4
    varargout{1}=PETpve;
end

return
%_______________________________________________________________________
function nits = get_nits(nits)
% Get number of nonlinear iterations
% FORMAT nits = get_nits(nits)

if prod(nits) > 0,
    tmp2 = [1 3 5 8 12 16];
    tmp = find(tmp2 == nits);
    if isempty(tmp) tmp = length(tmp2); end;
    nits = spm_input(['# Nonlinear Iterations?'],'+1','m',...
        ['1  nonlinear iteration |3  nonlinear iterations'...
        '|5  nonlinear iterations|8  nonlinear iterations'...
        '|12 nonlinear iterations|16 nonlinear iterations'],tmp2, tmp);
else, nits = 0; end;
return;
%_______________________________________________________________________

%_______________________________________________________________________
function cutoff = get_cutoff(cutoff)
% Get cutoff frequency of DCT bases
% FORMAT cutoff = get_cutoof(cutoff)

tmp2 = [15 20 25 30 35 40 45 50 60 70 80 Inf];
tmp = find(tmp2 == cutoff);
if isempty(tmp) tmp = length(tmp2); end;
cutoff = spm_input('Cutoff spatial normalization','+1','m',...
    ['15mm cutoff|20mm cutoff|25mm cutoff|30mm cutoff|'...
    '35mm cutoff|40mm cutoff|45mm cutoff|50mm cutoff|'...
    '60mm cutoff|70mm cutoff|80mm cutoff|Affine only'], tmp2, tmp);
return;
%_______________________________________________________________________

%_______________________________________________________________________
function reg = get_reg(reg)
% Get amount of regularisation
% FORMAT reg = get_reg(reg)

tmp2 = [100 10 1 0.1 0.01];
tmp = find(tmp2 == reg);
if isempty(tmp) tmp = length(tmp2); end;
reg = spm_input('Nonlinear Regularization','+1','m',...
    ['Extremely Heavy regularization (100)|Heavy regularization (10)|'...
    'Medium regularization (1)|Light regularization (0.1)|'...
    'Very Light regularization (0.01)'], tmp2, tmp);
return;
%_______________________________________________________________________
%_______________________________________________________________________

function Vo = ApplyBrainMask(V,Mask)
% P =spm_get(Inf,'IMAGE',['Select images to Maskify']);
% V   = spm_vol(P);%spm_get(Inf,'*.IMAGE',['Select ',GroupLabel, ' images to normalize & segment']));
f = 'i1.*i2';
for i=1:length(V)
    %V=spm_vol(P(i));
    [pth,nam,ext] = fileparts(V(i).fname);
    Q = ['b' nam];

    %-Work out filename for output image
    %------------------------------------------------------------------
    Qdir = spm_str_manip(Q,'Hv');
    Qfil = [spm_str_manip(Q,'stv'),'.img'];
    if ~exist(Qdir,'dir')
        warning('Invalid directory: writing to current directory')
        Qdir = '.';
    end
    Q = spm_get('CPath',Qfil,Qdir);
    Vo = struct(	'fname',	Q,...
        'dim',		[V(i).dim(1:end)],...
        'mat',		V(i).mat,...
        'descrip',	'spm - algebra');
    tmp(1)=V;
    tmp(2)=Mask;
    Vo = spm_imcalc(tmp,Vo,f);
end

%_______________________________________________________________________
function weight = get_weight(weight,type)
% Get cutoff frequency of DCT bases
% FORMAT cutoff = get_cutoof(cutoff)
if (strcmp(type,'NRM template') | strcmp(type,'SEG template'))
    tmp2 = [0 1];
    choices = 'no mask|generic mask';
elseif (strcmp(type,'SEG source') | strcmp(type,'NRM source'))
    tmp2 = [0 1];
    choices = 'no mask|subject mask';

else
    error('incorrect type')
end
tmp = find(tmp2 == weight);
if isempty(tmp) tmp = length(tmp2); end;
weight = spm_input([type ' Weighting Mask'],'+1','m',...
    choices, tmp2, tmp);

return;
%_______________________________________________________________________

%
% %_______________________________________________________________________
% function mask = get_mask(varargin)
% % Get cutoff frequency of DCT bases
% % FORMAT cutoff = get_cutoof(cutoff)
% if (nargin > 0), pth = varargin{1}, end
% if (nargin > 1), patid = varargin{2}, end
%
%         mask = fullfile(pth,'edits','bin', ['MSK_',patid,'.img']);
% return;
% %_______________________________________________________________________


%_______________________________________________________________________
function ApplyMask = get_ApplyMask(ApplyMask)
%
% Get cutoff frequency of DCT bases
% FORMAT cutoff = get_cutoof(cutoff)
tmp2 = [0 1];
choices = 'no, don''t apply mask|yes, apply mask';
tmp = find(tmp2 == ApplyMask);
if isempty(tmp) tmp = length(tmp2); end;
ApplyMask = spm_input(['Apply binary mask after segmentation?'],'+1','m',...
    choices, tmp2, tmp);

return;
%___________________________________________________________________
% function %%%%mls_cleanup(imagename)
% for i=1:length(imagename)
%     [pth,nam,ext,ver]=fileparts(imagename{i});
%     tmp1=fullfile(pth, [nam '.*' ver]);
%     disp(['clearing temporary file: ',imagename{i}])
%     delete(tmp1);
% end
