Unverified Commit ffee8fb5 authored by jayd1860's avatar jayd1860 Committed by GitHub
Browse files

Merge pull request #25 from BUNPC/multi-params 

hmrR_OD2Conc functions consolidated; hmrR_MotionArtifact & hmrR_MotionArtifactByChannel use mlActAuto; hmrR_tCCA input format fix
parents f0b7e845 2eee17d4

DataTree/ProcStream/ProcStreamClass.m

+1 −1
Original line number Diff line number Diff line
@@ -1181,7 +1181,7 @@ classdef ProcStreamClass < handle 
                    obj.reg.funcReg(iR).GetUsageStrDecorated(['hmrR_Intensity2OD',suffix]); ...

                    obj.reg.funcReg(iR).GetUsageStrDecorated(['hmrR_BandpassFilt',suffix],'aux'); ...

                    obj.reg.funcReg(iR).GetUsageStrDecorated(['hmrR_BandpassFilt',suffix],'dod'); ...

                    obj.reg.funcReg(iR).GetUsageStrDecorated(['hmrR_OD2Conc_new',suffix]); ...

                    obj.reg.funcReg(iR).GetUsageStrDecorated(['hmrR_OD2Conc',suffix]); ...

                    obj.reg.funcReg(iR).GetUsageStrDecorated(['hmrR_BlockAvg',suffix],'dcAvg'); ...

                };

                k=[]; kk=1;

FuncRegistry/UserFunctions/Archive/hmrR_OD2Conc_Nirs.m

+32 −14
Original line number Diff line number Diff line 
% SYNTAX:

% dc = hmrR_OD2Conc_Nirs( dod, SD, ppf )

% dc = hmrR_OD2Conc_new_Nirs( dod, SD, ppf )

%

% UI NAME:

% OD_to_Conc
@@ -14,36 +14,50 @@ 
%      of data, then this is a vector of 2 elements.  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: the concentration data (#time points x 3 x #SD pairs

%     3 concentrations are returned (HbO, HbR, HbT)

%

% USAGE OPTIONS:

% Delta_OD_to_Conc: dc = hmrR_OD2Conc_Nirs( dod, SD, ppf )

% Delta_OD_to_Conc: dc = hmrR_OD2Conc_new_Nirs( dod, SD, ppf )

%

% PARAMETERS:

% ppf: [1.0, 1.0]

% ppf: [1.0, 1.0, 1.0]

%

function dc = hmrR_OD2Conc_Nirs( dod, SD, ppf )

function dc = hmrR_OD2Conc_new_Nirs( dod, SD, ppf )



dc = [];

nWav = length(SD.Lambda);

ml = SD.MeasList;



if length(ppf)~=nWav

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

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

    return

if length(ppf) < nWav

    errordlg('WARNING: Data contains more than 3 wavelengths. Using PPF value of 1 for all wavelengths.');

    ppf = ones(1, nWav);

elseif length(ppf) > nWav

    d = length(ppf)-nWav;

    ppf(end-d+1:end) = [];

end



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

    ppf = ones(size(ppf));

end



nTpts = size(dod,1);



e = GetExtinctions( SD.Lambda );

if ~isfield(SD,'SpatialUnit')

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

elseif strcmpi(SD.SpatialUnit,'mm')

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

elseif strcmpi(SD.SpatialUnit,'cm')

    e = e(:,1:2) ;

end

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



lst = find( ml(:,4)==1 );
@@ -51,6 +65,10 @@ for idx=1:length(lst) 
    idx1 = lst(idx);

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

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

    if ppf(1)~=1

        dc(:,:,idx) = ( einv * (dod(:,[idx1 idx2'])./(ones(nTpts,1)*rho*ppf))' )';

    else

        dc(:,:,idx) = ( einv * (dod(:,[idx1 idx2'])./(ones(nTpts,1)))' )';

    end

end

dc(:,3,:) = dc(:,1,:) + dc(:,2,:);

FuncRegistry/UserFunctions/Archive/hmrR_OD2Conc_new_Nirs.m

deleted100644 → 0
+0 −71
Original line number Diff line number Diff line 
% SYNTAX:

% dc = hmrR_OD2Conc_new_Nirs( dod, SD, ppf )

%

% UI NAME:

% OD_to_Conc

%

% DESCRIPTION:

% Convert OD to concentrations

%

% INPUTS:

% dod: the change in OD (#time points x #channels)

% SD:  the SD structure

% 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 

%      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: the concentration data (#time points x 3 x #SD pairs

%     3 concentrations are returned (HbO, HbR, HbT)

%

% USAGE OPTIONS:

% Delta_OD_to_Conc: dc = hmrR_OD2Conc_new_Nirs( dod, SD, ppf )

%

% PARAMETERS:

% ppf: [1.0, 1.0]

%

function dc = hmrR_OD2Conc_new_Nirs( dod, SD, ppf )



nWav = length(SD.Lambda);

ml = SD.MeasList;



if length(ppf)~=nWav

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

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

    return

end



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

    ppf = ones(size(ppf));

end



nTpts = size(dod,1);



e = GetExtinctions( SD.Lambda );

if ~isfield(SD,'SpatialUnit')

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

elseif strcmpi(SD.SpatialUnit,'mm')

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

elseif strcmpi(SD.SpatialUnit,'cm')

    e = e(:,1:2) ;

end

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



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

for idx=1:length(lst)

    idx1 = lst(idx);

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

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

    if ppf(1)~=1

        dc(:,:,idx) = ( einv * (dod(:,[idx1 idx2'])./(ones(nTpts,1)*rho*ppf))' )';

    else

        dc(:,:,idx) = ( einv * (dod(:,[idx1 idx2'])./(ones(nTpts,1)))' )';

    end

end

dc(:,3,:) = dc(:,1,:) + dc(:,2,:);

FuncRegistry/UserFunctions/hmrR_MotionArtifact.m

+12 −4
Original line number Diff line number Diff line 
% SYNTAX:

% tInc = hmrR_MotionArtifact(data, probe, mlActMan, tIncMan, tMotion, tMask, STDEVthresh, AMPthresh)

% tInc = hmrR_MotionArtifact(data, probe, mlActMan, mlActAuto, tIncMan, tMotion, tMask, STDEVthresh, AMPthresh)

%

% UI NAME:

% Motion_Artifacts
@@ -15,6 +15,8 @@ 
% probe:   SNIRF data structure probe, containing probe source/detector geometry

% mlActMan: Cell array of vectors, one for each time base in data, specifying 

%        active/inactive channels with 1 meaning active, 0 meaning inactive

% mlActAuto: Cell array of vectors, one for each time base in data, specifying 

%        active/inactive channels with 1 meaning active, 0 meaning inactive

% tIncMan: Cell array of vectors corresponding to the number of time bases in data. 

%          tIncMan has been manually excluded. 0-excluded. 1-included. Vector same length as d.

% tMotion: Check for signal change indicative of a motion artifact over
@@ -42,7 +44,7 @@ 
%       with 1's indicating data included and 0's indicate motion artifact

%

% USAGE OPTIONS:

% Motion_Artifacts:  tIncAuto = hmrR_MotionArtifact(dod, probe, mlActMan, tIncMan, tMotion, tMask, STDEVthresh, AMPthresh)

% Motion_Artifacts:  tIncAuto = hmrR_MotionArtifact(dod, probe, mlActMan, mlActAuto, tIncMan, tMotion, tMask, STDEVthresh, AMPthresh)

%

% PARAMETERS:

% tMotion: 0.5
@@ -60,7 +62,7 @@ 
% JDUBB 3/18/2019 Adapted to SNIRF format

%

%

function tInc = hmrR_MotionArtifact(data, probe, mlActMan, tIncMan, tMotion, tMask, std_thresh, amp_thresh)

function tInc = hmrR_MotionArtifact(data, probe, mlActMan, mlActAuto, tIncMan, tMotion, tMask, std_thresh, amp_thresh)



% Init output 

tInc = cell(length(data),1);
@@ -81,6 +83,9 @@ end 
if isempty(mlActMan)

    mlActMan = cell(length(data),1);

end

if isempty(mlActAuto)

    mlActAuto = cell(length(data),1);

end



for iBlk=1:length(data)
@@ -100,7 +105,10 @@ for iBlk=1:length(data) 
    if isempty(mlActMan{iBlk})

        mlActMan{iBlk} = ones(size(MeasList,1),1);

    end

    MeasListAct = mlActMan{iBlk};        

    if isempty(mlActAuto{iBlk})

        mlActAuto{iBlk} = ones(size(MeasList,1),1);

    end

    MeasListAct = mlActMan{iBlk} & mlActAuto{iBlk};        

    

    % Calculate the diff of d to to set the threshold if ncssesary

    diff_d = diff(d);

FuncRegistry/UserFunctions/hmrR_MotionArtifactByChannel.m

+12 −4
Original line number Diff line number Diff line 
% [tInc,tIncCh] = hmrR_MotionArtifactByChannel(data, probe, mlActMan, tIncMan, tMotion, tMask, STDEVthresh, AMPthresh)

% [tInc,tIncCh] = hmrR_MotionArtifactByChannel(data, probe, mlActMan, mlActAuto, tIncMan, tMotion, tMask, STDEVthresh, AMPthresh)

%

% UI NAME:   

% Motion_Artifacts_By_Channel
@@ -19,6 +19,8 @@ 
%          tIncMan has been manually excluded. 0-excluded. 1-included. Vector same length as d.

% mlActMan: Cell array of vectors, one for each time base in data, specifying 

%        active/inactive channels with 1 meaning active, 0 meaning inactive

% mlActAuto: Cell array of vectors, one for each time base in data, specifying 

%        active/inactive channels with 1 meaning active, 0 meaning inactive

% tMotion: Check for signal change indicative of a motion artifact over

%     time range tMotion. Units of seconds.

% tMask: Mark data over +/- tMask seconds around the identified motion
@@ -40,7 +42,7 @@ 
%       channel basis

%

% USAGE OPTIONS:

% Motion_Artifacts_By_Channel:  [tIncAuto, tIncAutoCh] = hmrR_MotionArtifactByChannel(dod, probe, mlActMan, tIncMan, tMotion, tMask, STDEVthresh, AMPthresh)

% Motion_Artifacts_By_Channel:  [tIncAuto, tIncAutoCh] = hmrR_MotionArtifactByChannel(dod, probe, mlActMan, mlActAuto, tIncMan, tMotion, tMask, STDEVthresh, AMPthresh)

%

% PARAMETERS:

% tMotion: 0.5
@@ -59,7 +61,7 @@ 
% TO DO:

% Consider tIncMan



function [tInc, tIncCh] = hmrR_MotionArtifactByChannel(data, probe, mlActMan, tIncMan, tMotion, tMask, std_thresh, amp_thresh)

function [tInc, tIncCh] = hmrR_MotionArtifactByChannel(data, probe, mlActMan, mlActAuto, tIncMan, tMotion, tMask, std_thresh, amp_thresh)



tInc   = cell(length(data), 1);

tIncCh = cell(length(data), 1);
@@ -80,6 +82,9 @@ end 
if isempty(mlActMan)

    mlActMan = cell(length(data),1);

end

if isempty(mlActAuto)

    mlActAuto = cell(length(data),1);

end



for iBlk=1:length(data)
@@ -102,7 +107,10 @@ for iBlk=1:length(data) 
    if isempty(mlActMan{iBlk})

        mlActMan{iBlk} = ones(size(MeasList,1),1);

    end

    MeasListAct = mlActMan{iBlk};

    if isempty(mlActAuto{iBlk})

        mlActAuto{iBlk} = ones(size(MeasList,1),1);

    end

    MeasListAct = mlActMan{iBlk} & mlActAuto{iBlk};      

        

    % Calculate the diff of d to to set the threshold if ncssesary

    diff_d = diff(d);

CRA Git | Maintained and supported by SUSTech CRA and CCSE