Merge pull request #1 from BUNPC/master

                ml(ii,:) = [obj.measurementList(ii).GetSourceIndex(), obj.measurementList(ii).GetDetectorIndex(), 1, obj.measurementList(ii).GetWavelengthIndex()];

            end

            % Remove any unused allocated rows that were pre-allocated

            % Remove unused rows that were pre-allocated

            ml(ii+1:end,:) = [];

            % Sort according to wavelength

            ml = sortrows(ml,4);

        end

            obj.input.misc = [];

            close(hwait);

        end

            %            argOut: '[dcAvg,dcAvgStd,tHRF,nTrials]'

            %             argIn: '(dcAvgRuns,dcAvgStdRuns,dcSum2Runs,tHRFRuns,mlActRuns,nTrialsRuns'

            %           paramIn: [0x0 ParamClass]

            %              help: '  Calculate the block average for all subjects, for all common stimuli'

            %              help: '  Calculate the block average for all subjects, for all common stimuli�'

            %

            err = -1;

            if ~exist('fname', 'var')

            %    Here's the output:

            %

            %     G = {

            %          '@ hmrG_SubjAvg [dcAvg,dcAvgStd,nTrials,grpAvgPass] (dcAvgSubjs,dcAvgStdSubjs,SDSubjs,nTrialsSubjs tRange %0.1f'

            %          '@ hmrG_SubjAvg [dcAvg,dcAvgStd,nTrials,grpAvgPass] (dcAvgSubjs,dcAvgStdSubjs,SDSubjs,nTrialsSubjs tRange %0.1f�'

            %         }

            %     S = {

            %          '@ hmrS_RunAvg [dcAvg,dcAvgStd,nTrials] (dcAvgRuns,dcAvgStdRuns,dcSum2Runs,mlActRuns,nTrialsRuns'

            %         '@ hmrR_MotionArtifact tIncAuto (dod,t,SD,tIncMan tMotion %0.1f 0.5 tMask %0.1f 1.0 STDEVthresh %0.1f 50.0 AMPthresh %0.1f 5.0'

            %         '@ hmrR_BandpassFilt dod (dod,t hpf %0.3f 0.010 lpf %0.3f 0.500'

            %         '@ hmrR_OD2Conc dc (dod,SD ppf %0.1f_%0.1f 6.0_6.0'

            %         '@ hmrR_DeconvHRF_DriftSS [dcAvg,dcAvgstd,tHRF,nTrials,ynew,yresid,ysum2,beta,R] (dc,s,t,SD,aux,tIncAuto trange %0.1f_%0.1f -2.0_20.0 glmSolv'

            %         '@ hmrR_DeconvHRF_DriftSS [dcAvg,dcAvgstd,tHRF,nTrials,ynew,yresid,ysum2,beta,R] (dc,s,t,SD,aux,tIncAuto trange %0.1f_%0.1f -2.0_20.0 glmSolv�'

            %         }

            % 

            if ~exist('mode','var') || isempty(mode) || ~ischar(mode)

            %         '@ hmrR_MotionArtifact tIncAuto (dod,t,SD,tIncMan tMotion %0.1f 0.5 tMask %0.1f 1.0 STDEVthresh %0.1f 50.0 AMPthresh %0.1f 5.0'

            %         '@ hmrR_BandpassFilt dod (dod,t hpf %0.3f 0.010 lpf %0.3f 0.500'

            %         '@ hmrR_OD2Conc dc (dod,SD ppf %0.1f_%0.1f 6.0_6.0'

            %         '@ hmrR_DeconvHRF_DriftSS [dcAvg,dcAvgstd,tHRF,nTrials,ynew,yresid,ysum2,beta,R] (dc,s,t,SD,aux,tIncAuto trange %0.1f_%0.1f -2.0_20.0 glmSolv'

            %         '@ hmrR_DeconvHRF_DriftSS [dcAvg,dcAvgstd,tHRF,nTrials,ynew,yresid,ysum2,beta,R] (dc,s,t,SD,aux,tIncAuto trange %0.1f_%0.1f -2.0_20.0 glmSolv�'

            %         }

            %

            section = cell(length(obj.fcalls), 1);

            %            argOut: 'tIncAuto'

            %             argIn: '(dod,t,SD,tIncMan'

            %           paramIn: [1x4 ParamClass]

            %              help: '  Excludes stims that fall within the time points identified as '

            %              help: '  Excludes stims that fall within the time points identified as �'

            %

            %    p.fcalls(3)

            %     

            %            argOut: '[s,tRangeStimReject]'

            %             argIn: '(t,s,tIncAuto,tIncMan'

            %           paramIn: [1x1 ParamClass]

            %              help: '  Excludes stims that fall within the time points identified as '

            %              help: '  Excludes stims that fall within the time points identified as �'

            %

            err=-1;

            %        argOut: 'dod'

            %         argIn: '(dod,t'

            %       paramIn: [1x2 ParamClass]

            %          help: '  Perform a bandpass filter'

            %          help: '  Perform a bandpass filter�'

            %

            %    p.fcalls(2)

            %

            %        argOut: 'dc'

            %         argIn: '(dod,SD'

            %       paramIn: [1x1 ParamClass]

            %          help: '  Convert OD to concentrations'

            %          help: '  Convert OD to concentrations�'

            % 

            if nargin<2

                return

% OD_to_Conc

%

% DESCRIPTION:

% Convert OD to concentrations

% Convert OD to concentrations. For use with 2-wavelength data only: use

% hmrR_OD2Conc_3 or hmrR_OD2Conc_multi for 3 or more wavelengths.

%

% INPUTS:

% dod: SNIRF.data container with the Change in OD tim course 

% probe: SNIRF.probe container with the source/detector geometry

% ppf: Partial path length factors for each wavelength. If there are 2 wavelengths 

%      of data, then this is a vector of 2 elements.  Typical value is ~6 for each 

% ppf: Partial path length factors for each wavelength. This is a vector of 

%      of 2 elements, one for each wavelength.  Typical value is ~6 for each 

%      wavelength if the absorption change is uniform over the volume of tissue measured. 

%      To approximate the partial volume effect of a small localized absorption change 

%      within an adult human head, this value could be as small as 0.1. It is recommended 

%      to use default values of 1 1 which will result in concentration units of 

%      molar ppf such that the user can then divide by an estimated ppf at any future 

%      point to estimate what the molar concentration change would be.%

%      within an adult human head, this value could be as small as 0.1. Convention is 

%      becoming to set ppf=1 and to not divide by the source-detector separation such that 

%      the resultant "concentration" is in units of Molar mm (or Molar cm if those are the 

%      spatial units). This is becoming wide spread in the literature but there is no 

%      fixed citation. Use a value of 1 to choose this option.

%

% OUTPUTS:

% dc: SNIRF.data container with the concentration data 

% SYNTAX:

% dc = hmrR_OD2Conc_new( dod, probe, ppf )

%

% UI NAME:

% OD_to_Conc

%

% DESCRIPTION:

% Convert OD to concentrations. For use with 3-wavelength data only.

%

% INPUTS:

% dod: SNIRF.data container with the Change in OD tim course 

% probe: SNIRF.probe container with the source/detector geometry

% ppf: Partial path length factors for each wavelength. This is a vector of  

%      3 elements, one for each wavelength.  Typical value is ~6 for each 

%      wavelength if the absorption change is uniform over the volume of tissue measured. 

%      To approximate the partial volume effect of a small localized absorption change 

%      within an adult human head, this value could be as small as 0.1. Convention is 

%      becoming to set ppf=1 and to not divide by the source-detector separation such that 

%      the resultant "concentration" is in units of Molar mm (or Molar cm if those are the 

%      spatial units). This is becoming wide spread in the literature but there is no 

%      fixed citation. Use a value of 1 to choose this option.

%

% OUTPUTS:

% dc: SNIRF.data container with the concentration data 

%

% USAGE OPTIONS:

% Delta_OD_to_Conc: dc = hmrR_OD2Conc_new( dod, probe, ppf )

%

% PARAMETERS:

% ppf: [1.0, 1.0, 1.0]

%

function dc = hmrR_OD2Conc_new( dod, probe, ppf )

dc = DataClass().empty();

for ii=1:length(dod)

    dc(ii) = DataClass();

    Lambda = probe.GetWls();

    SrcPos = probe.GetSrcPos();

    DetPos = probe.GetDetPos();

    nWav   = length(Lambda);

    ml     = dod(ii).GetMeasList();

    y      = dod(ii).GetDataTimeSeries();

    if length(ppf)~=nWav

        errordlg('The length of PPF must match the number of wavelengths in SD.Lambda');

        dc = zeros(size(y,1),3,length(find(ml(:,4)==1)));

        return

    end

    if ~isempty(find(ppf==1))

        ppf = ones(size(ppf));

    end

    nTpts = size(y,1);

    e = GetExtinctions(Lambda);

    e = e(:,1:2) / 10; % convert from /cm to /mm

    einv = inv( e'*e )*e';

    lst = find( ml(:,4)==1 );

    y2 = zeros(nTpts, 3*length(lst));

    for idx=1:length(lst)

        k = 3*(idx-1)+1;

        idx1 = lst(idx);

        idx2 = find( ml(:,4)>1 & ml(:,1)==ml(idx1,1) & ml(:,2)==ml(idx1,2) );

        rho = norm(SrcPos(ml(idx1,1),:)-DetPos(ml(idx1,2),:));

        if ppf(1)~=1

            y2(:,k:k+1) = ( einv * (y(:,[idx1 idx2'])./(ones(nTpts,1)*rho*ppf))' )';

        else

            y2(:,k:k+1) = ( einv * (y(:,[idx1 idx2'])./(ones(nTpts,nWav)))' )';

        end

        y2(:,k+2) = y2(:,k) + y2(:,k+1);

        dc(ii).AddChannelHbO(ml(idx1,1), ml(idx1,2));

        dc(ii).AddChannelHbR(ml(idx1,1), ml(idx1,2));

        dc(ii).AddChannelHbT(ml(idx1,1), ml(idx1,2));

    end   

    dc(ii).SetDataTimeSeries(y2);

    dc(ii).SetTime(dod(ii).GetTime());

end