From c90b4ad818fe1f0055de9961b996f9b70834c38a Mon Sep 17 00:00:00 2001 From: kk1995 Date: Wed, 4 Oct 2023 17:18:50 -0400 Subject: [PATCH 01/20] Update setpaths.m Updated to be more generalizable to different path formats. For example some paths have '\\' at the beginning and previous setpaths did not support such paths. This one should support all paths supported by matlab via fullfile() function. --- setpaths.m | 60 ++++++++++++++++++++++++++---------------------------- 1 file changed, 29 insertions(+), 31 deletions(-) diff --git a/setpaths.m b/setpaths.m index 8e62f7e..a491bca 100644 --- a/setpaths.m +++ b/setpaths.m @@ -46,7 +46,7 @@ function setpaths(addremove) appNameInclList = {'Homer3'}; exclSearchList = {'.git','.idea','Data','Docs','*_install','*.app','submodules'}; - appThis = filesepStandard_startup(pwd); + appThis = pwd; appThisPaths = findDotMFolders(appThis, exclSearchList); if addremove == 0 @@ -224,15 +224,14 @@ function removeSearchPaths(app) % ---------------------------------------------------- function d = addDependenciesSearchPaths() -if exist([pwd, '/Utils/submodules'],'dir') - addpath([pwd, '/Utils/submodules'],'-end'); +if exist(fullfile(pwd, 'Utils','submodules'),'dir') + addpath(fullfile(pwd, 'Utils','submodules'),'-end'); end d = dependencies(); for ii = 1:length(d) rootpath = findFolder(pwd, d{ii}); - rootpath(rootpath=='\') = '/'; - if ispathvalid_startup([rootpath, '/Shared'],'dir') - rootpath = [rootpath, '/Shared']; + if ispathvalid_startup(fullfile(rootpath, 'Shared'),'dir') + rootpath = fullfile(rootpath, 'Shared'); end if ~exist(rootpath,'dir') printMethod(sprintf('ERROR: Could not find required dependency %s\n', d{ii})); @@ -362,12 +361,12 @@ function printMethod(msg) function dirpath = findFolder(repo, dirname) dirpath = ''; if ~exist('repo','var') - repo = filesepStandard_startup(pwd); + repo = pwd; end dirpaths = findDotMFolders(repo, {'.git', '.idea'}); for ii = 1:length(dirpaths) - [~, f, e] = fileparts(dirpaths{ii}(1:end-1)); + [~, f, e] = fileparts(dirpaths{ii}); if strcmp(dirname, [f,e]) dirpath = dirpaths{ii}; break; @@ -399,9 +398,10 @@ function printMethod(msg) exclList = {exclList}; end -subdirFullpath = filesepStandard_startup(subdir,'full'); +% subdirFullpath = filesepStandard_startup(subdir,'full'); +subdirFullpath = fullfile(subdir); -if ~ispathvalid_startup(subdirFullpath, 'dir') +if exist(subdirFullpath) ~= 7 logger.Write('Warning: folder %s doesn''t exist\n', subdirFullpath); return; end @@ -411,13 +411,14 @@ function printMethod(msg) return; end -dirs = dir([subdirFullpath, '*']); +dirs = dir(fullfile(subdirFullpath, '*')); +dirs = dirs(3:end); % ignore . and .. if isempty(dirs) return; end if isdotmfolder(subdirFullpath) - dotmfolders = {filesepStandard_startup(subdirFullpath, 'nameonly')}; + dotmfolders = {subdirFullpath}; end for ii = 1:length(dirs) @@ -427,7 +428,7 @@ function printMethod(msg) if dirs(ii).name(1) == '.' continue; end - dotmfolders = [dotmfolders; findDotMFolders([subdirFullpath, dirs(ii).name], exclList)]; %#ok + dotmfolders = [dotmfolders; findDotMFolders(fullfile(subdirFullpath, dirs(ii).name), exclList)]; %#ok end @@ -441,7 +442,7 @@ function printMethod(msg) if ~ispathvalid_startup(folder, 'dir') return end -if isempty(dir([folder,'/*.m'])) +if isempty(dir(fullfile(folder,'*.m'))) % Exceptions to rule that 'dotm' folder must have at least one '.m' file: % it is a an executable folder (i.e. '/bin') if ~isempty(strfind(folder, '/bin/')) %#ok<*STREMP> @@ -450,15 +451,16 @@ function printMethod(msg) end return; else - rootdir = which('findDotMFolders'); - rootdir = fileparts(rootdir); - rootdir = pathsubtract(rootdir, 'Utils/submodules','nochange'); - p = pathsubtract(folder, rootdir); - if length(find(p=='/')) > MAXPATHLENGTH - return - end + b = true; +% rootdir = which('findDotMFolders'); +% rootdir = fileparts(rootdir); +% rootdir = pathsubtract(rootdir, fullfile('Utils','submodules'),'nochange'); +% p = pathsubtract(folder, rootdir); +% if length(find(p==filesep)) > MAXPATHLENGTH +% return +% end end -b = true; +% b = true; @@ -496,7 +498,7 @@ function printMethod(msg) % ------------------------------------------------------------------------- -function diff = pathsubtract(p2_0, p1_0, options) +function diff = pathsubtract(p2, p1, options) if ~exist('options','var') options = ''; end @@ -505,14 +507,6 @@ function printMethod(msg) else option = 'full'; end -p1 = filesepStandard_startup(p1_0, option); -p2 = filesepStandard_startup(p2_0, option); -if isempty(p1) - p1 = p1_0; -end -if isempty(p2) - p2 = p2_0; -end k = strfind(p2, p1); if ~isempty(k) && k(1)==1 diff = p2(k(1)+length(p1):end); @@ -522,6 +516,10 @@ function printMethod(msg) diff = ''; end +if isempty(p1) + diff = p2; +end + % -------------------------------------------------------------------- From 16c742ab537e783282d7878ba0442cbb92f7470c Mon Sep 17 00:00:00 2001 From: kk1995 Date: Wed, 4 Oct 2023 17:36:23 -0400 Subject: [PATCH 02/20] removed filesepStandard_startup removed filesepStandard_startup() function from working on other parts of setpaths as it is not general enough compared to fullfile() --- setpaths.m | 41 ++++++++++++++++++++++------------------- 1 file changed, 22 insertions(+), 19 deletions(-) diff --git a/setpaths.m b/setpaths.m index a491bca..07c8c25 100644 --- a/setpaths.m +++ b/setpaths.m @@ -74,7 +74,7 @@ function setpaths(addremove) end for jj = 1:length(foo) - p = filesepStandard_startup(fileparts(foo{jj})); + p = fileparts(foo{jj}); if pathscompare_startup(appThis, p) continue end @@ -88,11 +88,11 @@ function setpaths(addremove) foo = which([appNameInclList{ii}, '.m'],'-all'); for jj = 1:length(foo) if jj > 1 - p = filesepStandard_startup(fileparts(foo{jj})); + p = fileparts(foo{jj}); appExclList = [appExclList; p]; %#ok printMethod(sprintf('Exclude paths for %s\n', p)); else - p = filesepStandard_startup(fileparts(foo{jj})); + p = fileparts(foo{jj}); appInclList = [appInclList; p]; %#ok printMethod(sprintf('Include paths for %s\n', p)); end @@ -188,7 +188,7 @@ function addSearchPaths(appPaths) addpath(appPaths{kk}, '-end'); setpermissions(appPaths{kk}); end -printMethod(sprintf('ADDED search paths for app %s\n', appPaths{1})); +printMethod(sprintf('ADDED search paths for app %s\n', string(appPaths{1}))); @@ -212,16 +212,13 @@ function removeSearchPaths(app) if ~isempty(strfind(lower(p{kk}), 'matlab')) && ~isempty(strfind(p{kk}, r)) continue; end - if ~isempty(strfind(filesepStandard_startup(p{kk}), app)) + if ~isempty(strfind(p{kk}, app)) rmpath(p{kk}); end end close(h); printMethod(sprintf('REMOVED search paths for app %s\n', app)); - - - % ---------------------------------------------------- function d = addDependenciesSearchPaths() if exist(fullfile(pwd, 'Utils','submodules'),'dir') @@ -240,9 +237,6 @@ function removeSearchPaths(app) addSearchPaths(rootpath); end - - - % ----------------------------------------------------------------------------- function [C,k] = str2cell_startup(str, delimiters, options) @@ -398,7 +392,6 @@ function printMethod(msg) exclList = {exclList}; end -% subdirFullpath = filesepStandard_startup(subdir,'full'); subdirFullpath = fullfile(subdir); if exist(subdirFullpath) ~= 7 @@ -781,13 +774,23 @@ function printMethod(msg) currdir = pwd; -% If paths are files, compare just the file names, then the folders -fullpath1 = filesepStandard_startup(path1, 'full'); -fullpath2 = filesepStandard_startup(path2, 'full'); - - -% Compare folders -b = strcmpi(fullpath1, fullpath2); +% If paths are files, compare just the file names, then the folders +if ~isfolder(path1) + [~,file1] = fileparts(path1); +else + file1 = []; +end +if ~isfolder(path2) + [~,file2] = fileparts(path2); +else + file2 = []; +end +if isempty(file1) && isempty(file2) + % Compare folders + b = strcmpi(path1, path2); +else + b = strcmpi(file1, file2); +end cd(currdir); From bff1df12074bdb88ff82f3a78f7177d76df22fc6 Mon Sep 17 00:00:00 2001 From: sreekanth kura Date: Wed, 11 Oct 2023 10:41:37 -0400 Subject: [PATCH 03/20] Fix an issue while changing Atlas --- AtlasViewerGUI.m | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/AtlasViewerGUI.m b/AtlasViewerGUI.m index 9c68835..d796a04 100644 --- a/AtlasViewerGUI.m +++ b/AtlasViewerGUI.m @@ -57,8 +57,10 @@ function ParseArgs(args) atlasViewer.dirnameAtlas = getAtlasDir(args); end if length(args)>3 - atlasViewer.handles.dataTree = args{4}; - atlasViewer.dataTree = get(atlasViewer.handles.dataTree, 'userdata'); + if sum(isgraphics(args{4})) + atlasViewer.handles.dataTree = args{4}; + atlasViewer.dataTree = get(atlasViewer.handles.dataTree, 'userdata'); + end end % Change current folder to dirnameSubj and load data From 20a507ec5d55dff4360dd7b377ea2727d18aa002 Mon Sep 17 00:00:00 2001 From: sreekanth kura Date: Mon, 30 Oct 2023 13:47:23 -0400 Subject: [PATCH 04/20] Fix bug in tranforming brain vol to head space in import MRI anatomy utility --- Utils/fs2viewer/segment_head_vols.m | 21 ++++++++++++++++++--- 1 file changed, 18 insertions(+), 3 deletions(-) diff --git a/Utils/fs2viewer/segment_head_vols.m b/Utils/fs2viewer/segment_head_vols.m index b9d1ffa..95aca6a 100644 --- a/Utils/fs2viewer/segment_head_vols.m +++ b/Utils/fs2viewer/segment_head_vols.m @@ -116,12 +116,27 @@ end % Align layer volume with head volume -X = inv(layer.vox2ras) * head.vox2ras; -layer.vol = xform_apply_vol_smooth(layer.vol, X); +X = inv(inv(layer.vox2ras) * head.vox2ras); +% layer.vol = xform_apply_vol_smooth(layer.vol, X); + +tform = affinetform3d(X); +% layer.vol = imwarp(layer.vol,tform,'OutputView',imref3d(size(head.vol))); % Assign layer segmentation number s = max(hseg(:))+1; -layer_i = find(layer.vol ~= 0); +% layer_i = find(layer.vol ~= 0); + +ind = find(layer.vol ~=0); +[I1,I2,I3] = ind2sub(size(layer.vol),ind); +t_ind = round([I1 I2 I3 ones(size(I1))]*X'); +idx = find(t_ind(:,1) < 1 | t_ind(:,1) > size(head.vol,1)); +t_ind(idx,:) = []; +idx = find(t_ind(:,2) < 1 | t_ind(:,2) > size(head.vol,2)); +t_ind(idx,:) = []; +idx = find(t_ind(:,3) < 1 | t_ind(:,3) > size(head.vol,3)); +t_ind(idx,:) = []; +layer_i = sub2ind(size(head.vol),t_ind(:,1),t_ind(:,2),t_ind(:,3)); + hseg(layer_i) = s; tiss_type{s}=name; From 9bc565ba504ffd4e24675540de1eda20b723db66 Mon Sep 17 00:00:00 2001 From: jayd1860 <44243610+jayd1860@users.noreply.github.com> Date: Thu, 16 Nov 2023 13:55:01 -0500 Subject: [PATCH 05/20] v2.48.1, Fix image reconstruction generating bad image that might be caused by inconsistent meas list channel order. (#65) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit * v2.48.0 -- Fix image reconstruction generating bad image that might be caused by inconsistent meas list channel order. Found by Meryem. Issue looks like in hmrConc2OD.m -- Fix issue in probe registration where AV does not detect presence of spring registration when it calls compatibleProbes.m from getProbe.m at load time. -- Fix issue in probe loading when searching SD files for compatible probes and it is not detected because of faulty algorithm in compatibleProbes.m. -- Fix in ImageRecon error when setting timerange, such that it is more than or less that tHRF. Fix is to get tHRF when ImageRecon GUI launches and use ceiln and floorn for min and max tHRF to make sure it is not out of range. Added floorn.m and ceiln.m to fix this. -- Add 2nd wavelength to get_tiss_prop.m to be the default rather than one wavelength and set default tiss ptop values for 760 nm and 850 nm as published in Sensors. 13 February 2023: "Impact of Anatomical Variability on Sensitivity Profile in fNIRS–MRI Integration" -- Fix failure to load full head probe fullhead_56x144_v3.SD or ninjaNIRS2022_2023-09-19-18-26-16.snirf files because SnirfClass.GetAuxDataMatrix gets an error when resampling with Matlab's resample function which does not handle a huge difference between aux sample rate and the data sample rate. To fix added local resample function resamplee.m to avoid this issue. DataTree, v1.17.3 -- Fix failure to load full head probe fullhead_56x144_v3.SD or ninjaNIRS2022_2023-09-19-18-26-16.snirf files because SnirfClass.GetAuxDataMatrix gets an error when resampling with Matlab's resample function which does not handle a huge difference between aux sample rate and the data sample rate. To fix added local resample function resamplee.m to avoid this issue. Utils, 1.8.1 -- Fix in ImageRecon error when setting timerange, such that it is more than or less that tHRF. Fix is to get tHRF when ImageRecon GUI launches and use ceiln and floorn for min and max tHRF to make sure it is not out of range. Added floorn.m and ceiln.m to fix this. -- Fix failure to load full head probe fullhead_56x144_v3.SD or ninjaNIRS2022_2023-09-19-18-26-16.snirf files because SnirfClass.GetAuxDataMatrix gets an error when resampling with Matlab's resample function which does not handle a huge difference between aux sample rate and the data sample rate. To fix added local resample function resamplee.m to avoid this issue. * v2.48.1 DataTree, v1.17.4 -- Fix probe not displaying in Homer3 because probe.sourcePos2D and probe.detectorPos2D are all zeros (which signifies empty but is not detected as empty) instead of actually empty. Also try to make units for generated 2D probe same as 3D by calling new method ProbeClass.ScaleProjection2D() -- Fix AuxClass eq function to check for different dataTimeSeries sizes when checking if 2 AuxClass objects are equal --- AtlasViewerGUI.m | 52 ++----- ForwardModel/genSensitivityProfile.m | 31 ++-- ForwardModel/getFwmodel.m | 28 +--- ForwardModel/get_tiss_prop.m | 136 ++++++++++-------- Group/DataTree/AcquiredData/Snirf/AuxClass.m | 6 + .../DataTree/AcquiredData/Snirf/ProbeClass.m | 33 +++++ .../DataTree/AcquiredData/Snirf/SnirfClass.m | 6 +- Group/DataTree/Version.txt | 2 +- ImgRecon/ImageRecon.m | 32 +++-- ImgRecon/hmrConc2OD.m | 46 ++++-- Probe/compatibleProbes.m | 2 +- Probe/findUniqueProbeFiles.m | 5 +- Probe/initProbe.m | 1 - Probe/isPreRegisteredProbe.m | 4 - Probe/pullProbeToHeadsurf.m | 2 +- Utils/Shared/Version.txt | 2 +- Utils/Shared/ceiln.m | 8 ++ Utils/Shared/floorn.m | 9 ++ Utils/Shared/getAppDir.m | 4 +- Utils/Shared/pretty_print_struct.m | 8 +- Utils/Shared/resamplee.m | 10 ++ Version.txt | 2 +- 22 files changed, 242 insertions(+), 187 deletions(-) create mode 100644 Utils/Shared/ceiln.m create mode 100644 Utils/Shared/floorn.m create mode 100644 Utils/Shared/resamplee.m diff --git a/AtlasViewerGUI.m b/AtlasViewerGUI.m index 9c68835..c9d0dea 100644 --- a/AtlasViewerGUI.m +++ b/AtlasViewerGUI.m @@ -653,27 +653,10 @@ function pushbuttonRegisterProbeToSurface_Callback(hObject, eventdata, handles) % Finish registration if isPreRegisteredProbe(probe, refpts) - % check if probe is already on the surface - probe_on_surface = true; - if ~isempty(atlasViewer.probe.registration.al) - for u = 1:size(atlasViewer.probe.registration.al,1) - ii = ismember(atlasViewer.refpts.labels,atlasViewer.probe.registration.al{u,2}); - refpt_pos = atlasViewer.refpts.pos(ii,:); - opt_pos = atlasViewer.probe.optpos_reg(atlasViewer.probe.registration.al{u,1},:); - dist_refpt_to_opt = sqrt(sum(refpt_pos-opt_pos).^2); - if dist_refpt_to_opt > 1 - probe_on_surface = false; - break; - end - end - end - % Register probe by simply pulling (or pushing) optodes toward surface % toward (or away from) center of head. - if ~ probe_on_surface probe = pullProbeToHeadsurf(probe, headobj); probe.hOptodesIdx = 1; - end else @@ -1281,15 +1264,13 @@ function menuItemGenerateMCInput_Callback(hObject, eventdata, handles) imgrecon = atlasViewer.imgrecon; hbconc = atlasViewer.hbconc; probe = atlasViewer.probe; -headvol = atlasViewer.headvol; pialsurf = atlasViewer.pialsurf; -headsurf = atlasViewer.headsurf; dirnameSubj = atlasViewer.dirnameSubj; axesv = atlasViewer.axesv; try if isempty(eventdata) || strcmp(eventdata.EventName,'Action') - fwmodel = genSensitivityProfile(fwmodel,probe,headvol,pialsurf,headsurf,dirnameSubj); + fwmodel = genSensitivityProfile(fwmodel, probe, dirnameSubj); if isempty(fwmodel.Adot) return; end @@ -1512,7 +1493,7 @@ function menuItemFs2Viewer_Callback(hObject, eventdata, handles) % -------------------------------------------------------------------- -function menuItemShowRefpts_Callback(hObject, eventdata, handles) +function menuItemShowRefpts_Callback(hObject, ~, ~) global atlasViewer switch(get(hObject, 'tag')) @@ -1592,7 +1573,6 @@ function pushbuttonCopyFigure_Callback(~, ~, ~) axis equal axis vis3d set(gca, 'unit','normalized'); -p = get(gca, 'position'); set(gca, 'unit','normalized', 'position', [.20, .30, .40, .40]); % colormap is a propery of figure not axes. Since we don't want to @@ -1666,10 +1646,10 @@ function menuItemLighting_Callback(hObject, eventdata, handles) axesv = atlasViewer.axesv; -if strcmp(get(hObject,'checked'), 'on'); +if strcmp(get(hObject,'checked'), 'on') set(hObject,'checked', 'off'); val=0; -elseif strcmp(get(hObject,'checked'), 'off'); +elseif strcmp(get(hObject,'checked'), 'off') set(hObject,'checked', 'on'); val=1; end @@ -1686,12 +1666,12 @@ function menuItemLighting_Callback(hObject, eventdata, handles) % -------------------------------------------------------------------- -function menuItemFindRefpts_Callback(hObject, eventdata, handles) +function menuItemFindRefpts_Callback(~, ~, ~) FindRefptsGUI(); % -------------------------------------------------------------------- -function menuProbePlacementVariation_Callback(hObject, eventdata, handles) +function menuProbePlacementVariation_Callback(~, ~, ~) plotProbePlacementVariation(); @@ -1958,7 +1938,7 @@ function menuCalcOptProp_Callback(~, ~, ~) set(hLabelsSurf,'FaceVertexCData',faceVertexCData); faceVertexAlphaData(iFace(iFaceMin)) = ones(length(iFace(iFaceMin)),1); set(hLabelsSurf,'FaceVertexAlphaData',faceVertexAlphaData); - iFaces = [iFaces iFace(iFaceMin)]; + iFaces = [iFaces, iFace(iFaceMin)]; end if all(ishandles(hProjectionRays)) set(hProjectionRays,'color','k'); @@ -2245,13 +2225,8 @@ function menuItemGetSensitivityatMNICoordinates_Callback(~, ~, handles) % -------------------------------------------------------------------- -function checkbox_Display_MNI_Projection_Callback(hObject, eventdata, handles) +function checkbox_Display_MNI_Projection_Callback(hObject, ~, ~) global atlasViewer - -fwmodel = atlasViewer.fwmodel; -headvol = atlasViewer.headvol; -labelssurf = atlasViewer.labelssurf; -headvol = atlasViewer.headvol; headsurf = atlasViewer.headsurf; if get(hObject,'value')==1 % if checkbox is checked @@ -2268,12 +2243,12 @@ function checkbox_Display_MNI_Projection_Callback(hObject, eventdata, handles) foo = num2str(coordinate_mni); fooi = sprintf(foo(1,:)); if no_mni == 1 - defaultanswer = {[fooi]}; + defaultanswer = {fooi}; else for i = 2:no_mni foon = ['; ' sprintf(foo(i,:))]; defaultanswer = {[fooi foon]}; - fooi = [fooi foon]; + fooi = [fooi, foon]; end %defaultanswer = {[sprintf(foo(1,:)) ';' sprintf(foo(2,:))]} end else @@ -2300,7 +2275,6 @@ function checkbox_Display_MNI_Projection_Callback(hObject, eventdata, handles) h2 = findobj('Marker','o'); set(h2,'Visible','off'); - headvol = atlasViewer.headvol; vertices = atlasViewer.labelssurf.mesh.vertices; % Project MNI in MC space to head surface and pial surface @@ -2379,7 +2353,6 @@ function menuItemSaveAnatomy_Callback(~, ~, ~) pialsurf = atlasViewer.pialsurf; labelssurf = atlasViewer.labelssurf; refpts = atlasViewer.refpts; -dirnameSubj = atlasViewer.dirnameSubj; saveHeadvol(headvol); saveHeadsurf(headsurf, headvol.T_2mc); @@ -2390,12 +2363,11 @@ function menuItemSaveAnatomy_Callback(~, ~, ~) % -------------------------------------------------------------------- -function menuItemLoadPrecalculatedProfile_Callback(hObject, eventdata, handles) +function menuItemLoadPrecalculatedProfile_Callback(~, ~, handles) global atlasViewer fwmodel = atlasViewer.fwmodel; headvol = atlasViewer.headvol; -pialsurf = atlasViewer.pialsurf; probe = atlasViewer.probe; T_vol2mc = headvol.T_2mc; @@ -2482,7 +2454,7 @@ function popupmenuImageDisplay_Callback(hObject, ~, handles) % -------------------------------------------------------------------- -function popupmenuImageDisplay_CreateFcn(hObject, eventdata, handles) +function popupmenuImageDisplay_CreateFcn(hObject, ~, ~) global popupmenuorder; popupmenuorder = struct(... diff --git a/ForwardModel/genSensitivityProfile.m b/ForwardModel/genSensitivityProfile.m index 16f4a61..64ebba4 100644 --- a/ForwardModel/genSensitivityProfile.m +++ b/ForwardModel/genSensitivityProfile.m @@ -1,4 +1,4 @@ -function fwmodel = genSensitivityProfile(fwmodel,probe,headvol,pialsurf,headsurf,dirnameSubj) +function fwmodel = genSensitivityProfile(fwmodel, probe, dirnameSubj) % Compute volume sensitivity matrix from Monte Carlo simulations and then % projects volume sensitivity onto pial mesh to obtain surface sensitivity @@ -12,6 +12,9 @@ % usage: % call function Adot in subject directory or registered atlas % directory. +global logger + +logger = InitLogger(logger); fwmodel = resetSensitivity(fwmodel,probe,dirnameSubj); if dirnameSubj(end)~='/' && dirnameSubj(end)~='\' @@ -76,13 +79,6 @@ Adot = single(zeros(nMeas,nNode,nWav)); % Adot_scalp = single(zeros(nMeas,nNode_scalp,nWav)); - - % We don't want to reserve the whole matrix in memory (too large). - % rather append each channel's sensitivity to a file. - A = single(zeros(nx,ny,nz)); - As = single(zeros(nx,ny,nz)); - Ad = single(zeros(nx,ny,nz)); - if fwmodel.AdotVolFlag fprintf('Warning: option to generate Adot 3 pt file enabled - May run out of memory\n'); fid1 = fopen([dirnameOut 'AdotVol.3pt'],'wb'); @@ -108,11 +104,11 @@ % load 2pt for given measurement from mc2 mcextreme output iS = probe.ml(iM,1); fileS = sprintf('%sfw%d.s%d.%s', dirnameOut, iW, iS, mc_output_ext); - As = loadMCFuncPtr(fileS, [nx ny nz 1]); + As = loadMCFuncPtr(fileS, [nx, ny, nz, 1]); iD = probe.ml(iM,2); fileD = sprintf('%sfw%d.d%d.%s', dirnameOut, iW, iD, mc_output_ext); - Ad = loadMCFuncPtr(fileD, [nx ny nz 1]); + Ad = loadMCFuncPtr(fileD, [nx, ny, nz, 1]); if fwmodel.normalizeFluence @@ -134,7 +130,7 @@ if sum_p~=0 As(idx_p) = As(idx_p) * (1 - abs(sum_n)) / sum_p; else - disp(sprintf('No photons launched into tissue from Src %d',iS)) + logger.Write('No photons launched into tissue from Src %d\n',iS) As(idx_p) = 0; end @@ -156,7 +152,7 @@ if sum_p~=0 Ad(idx_p) = Ad(idx_p) * (1 - abs(sum_n)) / sum_p; else - fprintf('No photons launched into tissue form Det %d\n',iD) + logger.Write('No photons launched into tissue form Det %d\n',iD) Ad(idx_p) = 0; end @@ -186,12 +182,18 @@ normfactor_d = As(floor(detpos{1}),floor(detpos{2}),floor(detpos{3})); end + [~, fs] = fileparts(fileS); + [~, fd] = fileparts(fileD); + logger.Write('Ch %d. Loading fluence files [%s, %s]\n', iM, fs, fd); + logger.Write('Generating sensitivity for ch %d = [%d, %d, %d], (s%d pos = [%0.1f, %0.1f, %0.1f], d%d pos = [%0.1f, %0.1f, %0.1f])\n', ... + iM, probe.ml(iM,1), probe.ml(iM,2), iW, iS, srcpos{1}, srcpos{2}, srcpos{3}, iD, detpos{1}, detpos{2}, detpos{3}); + % Get Adot normfactor = (normfactor_s + normfactor_d)/2; if normfactor~=0 A = (As.*Ad)/normfactor; else - fprintf('No photons detected between Src %d and Det %d',iS,iD); + logger.Write('No photons detected between Src %d and Det %d\n',iS,iD); A = zeros(size(As)); end if fwmodel.AdotVolFlag @@ -203,6 +205,9 @@ % memory function only available on windows memory; end + + + logger.Write('\n'); % Extract sensitivity values for the nodes of the low res mesh % (Adot) from the sensitivity matrix volume (A) diff --git a/ForwardModel/getFwmodel.m b/ForwardModel/getFwmodel.m index 893eaf7..c9517c5 100644 --- a/ForwardModel/getFwmodel.m +++ b/ForwardModel/getFwmodel.m @@ -123,33 +123,7 @@ % Look at MC input files for forward model params only % if not using precalculated fluence files if isempty(fwmodel.fluenceProfFnames) - - % Wavelength 1 - if exist([dirnameOut, 'fw1.s1.inp'],'file') - config = read_tMCimg_inp([dirnameOut, 'fw1.s1.inp']); - fwmodel.nphotons = config.phot_num; - for ii=1:size(config.tiss_prop,1) - fwmodel.headvol.tiss_prop(ii).scattering(:,1) = config.tiss_prop(ii,1); - fwmodel.headvol.tiss_prop(ii).anisotropy(:,1) = config.tiss_prop(ii,2); - fwmodel.headvol.tiss_prop(ii).absorption(:,1) = config.tiss_prop(ii,3); - fwmodel.headvol.tiss_prop(ii).refraction(:,1) = config.tiss_prop(ii,4); - end - fwmodel.nWavelengths = 1; - end - - % Wavelength 2 - if exist([dirnameOut, 'fw2.s1.inp'],'file') - config = read_tMCimg_inp([dirnameOut, 'fw2.s1.inp']); - fwmodel.nphotons = config.phot_num; - for ii=1:size(config.tiss_prop,1) - fwmodel.headvol.tiss_prop(ii).scattering(:,2) = config.tiss_prop(ii,1); - fwmodel.headvol.tiss_prop(ii).anisotropy(:,2) = config.tiss_prop(ii,2); - fwmodel.headvol.tiss_prop(ii).absorption(:,2) = config.tiss_prop(ii,3); - fwmodel.headvol.tiss_prop(ii).refraction(:,2) = config.tiss_prop(ii,4); - end - fwmodel.nWavelengths = 2; - end - + fwmodel.nWavelengths = length(fwmodel.headvol.tiss_prop(1).scattering); end end diff --git a/ForwardModel/get_tiss_prop.m b/ForwardModel/get_tiss_prop.m index 244176c..fffea06 100644 --- a/ForwardModel/get_tiss_prop.m +++ b/ForwardModel/get_tiss_prop.m @@ -46,38 +46,48 @@ tiss_prop = struct([]); -% Default tissue property values -SCATTERING_SKIN_DEF_VAL = 0.6600; -SCATTERING_SKULL_DEF_VAL = 0.8600; -SCATTERING_DM_DEF_VAL = 0.6600; -SCATTERING_CSF_DEF_VAL = 0.0100; -SCATTERING_GM_DEF_VAL = 1.1000; -SCATTERING_WM_DEF_VAL = 1.1000; -SCATTERING_OTHER_DEF_VAL = 0.8600; +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Default tissue properties for wavelength 760nm and 860nm +% Reference: Sensors. 13 February 2023 +% Impact of Anatomical Variability on Sensitivity Pro?le in fNIRS�MRI Integration +% Augusto Bonilauri 1, * , Francesca Sangiuliano Intra 2, Francesca Baglio 3and Giuseppe Baselli 1 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -ABSORPTION_SKIN_DEF_VAL = 0.0191; -ABSORPTION_SKULL_DEF_VAL = 0.0136; -ABSORPTION_DM_DEF_VAL = 0.0191; -ABSORPTION_CSF_DEF_VAL = 0.0026; -ABSORPTION_GM_DEF_VAL = 0.0186; -ABSORPTION_WM_DEF_VAL = 0.0186; -ABSORPTION_OTHER_DEF_VAL = 0.0191; +% Scattering +SCATTERING_SKIN_DEF_VAL = [0.6727, 0.5818]; +SCATTERING_SKULL_DEF_VAL = [0.8545, 0.7636]; +SCATTERING_DM_DEF_VAL = [0.6600, 0.2727]; +SCATTERING_CSF_DEF_VAL = [0.2727, 0.2727]; +SCATTERING_GM_DEF_VAL = [0.7600, 0.6165]; +SCATTERING_WM_DEF_VAL = [1.0825, 0.9188]; +SCATTERING_OTHER_DEF_VAL = [0.8600, 0.6727]; -ANISOTROPY_SKIN_DEF_VAL = 0.0010; -ANISOTROPY_SKULL_DEF_VAL = 0.0010; -ANISOTROPY_DM_DEF_VAL = 0.0010; -ANISOTROPY_CSF_DEF_VAL = 0.0010; -ANISOTROPY_GM_DEF_VAL = 0.0010; -ANISOTROPY_WM_DEF_VAL = 0.0010; -ANISOTROPY_OTHER_DEF_VAL = 0.0010; +% Absorption +ABSORPTION_SKIN_DEF_VAL = [0.0170, 0.0019]; +ABSORPTION_SKULL_DEF_VAL = [0.0116, 0.0139]; +ABSORPTION_DM_DEF_VAL = [0.0191, 0.6727]; +ABSORPTION_CSF_DEF_VAL = [0.0040, 0.0040]; +ABSORPTION_GM_DEF_VAL = [0.0180, 0.0192]; +ABSORPTION_WM_DEF_VAL = [0.0167, 0.0208]; +ABSORPTION_OTHER_DEF_VAL = [0.0180, 0.6727]; -REFRACTION_SKIN_DEF_VAL = 1.0000; -REFRACTION_SKULL_DEF_VAL = 1.0000; -REFRACTION_DM_DEF_VAL = 1.0000; -REFRACTION_CSF_DEF_VAL = 1.0000; -REFRACTION_GM_DEF_VAL = 1.0000; -REFRACTION_WM_DEF_VAL = 1.0000; -REFRACTION_OTHER_DEF_VAL = 1.0000; +% Anisotropy +ANISOTROPY_SKIN_DEF_VAL = [0.0010, 0.0010]; +ANISOTROPY_SKULL_DEF_VAL = [0.0010, 0.0010]; +ANISOTROPY_DM_DEF_VAL = [0.0010, 0.0010]; +ANISOTROPY_CSF_DEF_VAL = [0.0010, 0.0010]; +ANISOTROPY_GM_DEF_VAL = [0.0010, 0.0010]; +ANISOTROPY_WM_DEF_VAL = [0.0010, 0.0010]; +ANISOTROPY_OTHER_DEF_VAL = [0.0010, 0.0010]; + +% Refraction +REFRACTION_SKIN_DEF_VAL = [1.0000, 1.0000]; +REFRACTION_SKULL_DEF_VAL = [1.0000, 1.0000]; +REFRACTION_DM_DEF_VAL = [1.0000, 1.0000]; +REFRACTION_CSF_DEF_VAL = [1.0000, 1.0000]; +REFRACTION_GM_DEF_VAL = [1.0000, 1.0000]; +REFRACTION_WM_DEF_VAL = [1.0000, 1.0000]; +REFRACTION_OTHER_DEF_VAL = [1.0000, 1.0000]; %%% Extract args @@ -143,7 +153,7 @@ %%% Parse tissue names and tissue property names and find their values -propval = []; +propval = {}; m=length(propnames); n=length(names); for i=1:m @@ -152,84 +162,84 @@ for j=1:n switch lower(names{j}) case {'skin', 'scalp'} - propval(j,i) = ANISOTROPY_SKIN_DEF_VAL; + propval{j,i} = ANISOTROPY_SKIN_DEF_VAL; case {'skull', 'bone'} - propval(j,i) = ANISOTROPY_SKULL_DEF_VAL; + propval{j,i} = ANISOTROPY_SKULL_DEF_VAL; case {'dura', 'dura mater', 'dm'} - propval(j,i) = ANISOTROPY_DM_DEF_VAL; + propval{j,i} = ANISOTROPY_DM_DEF_VAL; case {'csf', 'cerebral spinal fluid'} - propval(j,i) = ANISOTROPY_CSF_DEF_VAL; + propval{j,i} = ANISOTROPY_CSF_DEF_VAL; case {'gm', 'gray matter', 'brain'} - propval(j,i) = ANISOTROPY_GM_DEF_VAL; + propval{j,i} = ANISOTROPY_GM_DEF_VAL; case {'wm', 'white matter'} - propval(j,i) = ANISOTROPY_WM_DEF_VAL; + propval{j,i} = ANISOTROPY_WM_DEF_VAL; case 'other' - propval(j,i) = ANISOTROPY_OTHER_DEF_VAL; + propval{j,i} = ANISOTROPY_OTHER_DEF_VAL; otherwise - propval(j,i) = ANISOTROPY_OTHER_DEF_VAL; + propval{j,i} = ANISOTROPY_OTHER_DEF_VAL; end end case 'scattering' for j=1:n switch lower(names{j}) case {'skin', 'scalp'} - propval(j,i) = SCATTERING_SKIN_DEF_VAL; + propval{j,i} = SCATTERING_SKIN_DEF_VAL; case {'skull', 'bone'} - propval(j,i) = SCATTERING_SKULL_DEF_VAL; + propval{j,i} = SCATTERING_SKULL_DEF_VAL; case {'dura', 'dura mater', 'dm'} - propval(j,i) = SCATTERING_DM_DEF_VAL; + propval{j,i} = SCATTERING_DM_DEF_VAL; case {'csf', 'cerebral spinal fluid'} - propval(j,i) = SCATTERING_CSF_DEF_VAL; + propval{j,i} = SCATTERING_CSF_DEF_VAL; case {'gm', 'gray matter', 'brain'} - propval(j,i) = SCATTERING_GM_DEF_VAL; + propval{j,i} = SCATTERING_GM_DEF_VAL; case {'wm', 'white matter'} - propval(j,i) = SCATTERING_WM_DEF_VAL; + propval{j,i} = SCATTERING_WM_DEF_VAL; case 'other' - propval(j,i) = SCATTERING_OTHER_DEF_VAL; + propval{j,i} = SCATTERING_OTHER_DEF_VAL; otherwise - propval(j,i) = SCATTERING_OTHER_DEF_VAL; + propval{j,i} = SCATTERING_OTHER_DEF_VAL; end end case 'absorption' for j=1:n switch lower(names{j}) case {'skin', 'scalp'} - propval(j,i) = ABSORPTION_SKIN_DEF_VAL; + propval{j,i} = ABSORPTION_SKIN_DEF_VAL; case {'skull', 'bone'} - propval(j,i) = ABSORPTION_SKULL_DEF_VAL; + propval{j,i} = ABSORPTION_SKULL_DEF_VAL; case {'dura', 'dura mater', 'dm'} - propval(j,i) = ABSORPTION_DM_DEF_VAL; + propval{j,i} = ABSORPTION_DM_DEF_VAL; case {'csf', 'cerebral spinal fluid'} - propval(j,i) = ABSORPTION_CSF_DEF_VAL; + propval{j,i} = ABSORPTION_CSF_DEF_VAL; case {'gm', 'gray matter', 'brain'} - propval(j,i) = ABSORPTION_GM_DEF_VAL; + propval{j,i} = ABSORPTION_GM_DEF_VAL; case {'wm', 'white matter'} - propval(j,i) = ABSORPTION_WM_DEF_VAL; + propval{j,i} = ABSORPTION_WM_DEF_VAL; case 'other' - propval(j,i) = ABSORPTION_OTHER_DEF_VAL; + propval{j,i} = ABSORPTION_OTHER_DEF_VAL; otherwise - propval(j,i) = ABSORPTION_OTHER_DEF_VAL; + propval{j,i} = ABSORPTION_OTHER_DEF_VAL; end end case 'refraction' for j=1:n switch lower(names{j}) case {'skin', 'scalp'} - propval(j,i) = REFRACTION_SKIN_DEF_VAL; + propval{j,i} = REFRACTION_SKIN_DEF_VAL; case {'skull', 'bone'} - propval(j,i) = REFRACTION_SKULL_DEF_VAL; + propval{j,i} = REFRACTION_SKULL_DEF_VAL; case {'dura', 'dura mater', 'dm'} - propval(j,i) = REFRACTION_DM_DEF_VAL; + propval{j,i} = REFRACTION_DM_DEF_VAL; case {'csf', 'cerebral spinal fluid'} - propval(j,i) = REFRACTION_CSF_DEF_VAL; + propval{j,i} = REFRACTION_CSF_DEF_VAL; case {'gm', 'gray matter', 'brain'} - propval(j,i) = REFRACTION_GM_DEF_VAL; + propval{j,i} = REFRACTION_GM_DEF_VAL; case {'wm', 'white matter'} - propval(j,i) = REFRACTION_WM_DEF_VAL; + propval{j,i} = REFRACTION_WM_DEF_VAL; case 'other' - propval(j,i) = REFRACTION_OTHER_DEF_VAL; + propval{j,i} = REFRACTION_OTHER_DEF_VAL; otherwise - propval(j,i) = REFRACTION_OTHER_DEF_VAL; + propval{j,i} = REFRACTION_OTHER_DEF_VAL; end end end @@ -240,6 +250,6 @@ for j=1:length(names) tiss_prop(j).name = names{j}; for i=1:length(propnames) - eval(sprintf('tiss_prop(j).%s = propval(j,i);',propnames{i})); + eval(sprintf('tiss_prop(j).%s = propval{j,i};',propnames{i})); end end diff --git a/Group/DataTree/AcquiredData/Snirf/AuxClass.m b/Group/DataTree/AcquiredData/Snirf/AuxClass.m index 028eed4..95bbb24 100644 --- a/Group/DataTree/AcquiredData/Snirf/AuxClass.m +++ b/Group/DataTree/AcquiredData/Snirf/AuxClass.m @@ -196,6 +196,12 @@ function Copy(obj, obj2) if ~strcmp(obj.name, obj2.name) return; end + if ~all(size(obj.dataTimeSeries)==size(obj2.dataTimeSeries)) + return; + end + if ~all(size(obj.time)==size(obj2.time)) + return; + end if ~all(obj.dataTimeSeries(:)==obj2.dataTimeSeries(:)) return; end diff --git a/Group/DataTree/AcquiredData/Snirf/ProbeClass.m b/Group/DataTree/AcquiredData/Snirf/ProbeClass.m index 451ded9..2c8fa91 100644 --- a/Group/DataTree/AcquiredData/Snirf/ProbeClass.m +++ b/Group/DataTree/AcquiredData/Snirf/ProbeClass.m @@ -150,8 +150,40 @@ function BackwardCompatibility(obj) + % ------------------------------------------------------- + function ScaleProjection2D(obj) + if size(obj.sourcePos2D,1) + dm_src2d = distmatrix([obj.sourcePos2D, zeros(size(obj.sourcePos2D,1),1)]); + else + dm_src2d = distmatrix(obj.sourcePos2D, zeros); + end + dm_src3d = distmatrix(obj.sourcePos3D); + exp = round(log10(mean(dm_src3d(dm_src3d(:)>0)))) - round(log10(mean(dm_src2d(dm_src2d(:)>0)))); + scaleFactor = 10^exp; + obj.sourcePos2D = obj.sourcePos2D*scaleFactor; + + if size(obj.detectorPos2D,1) + dm_det2d = distmatrix([obj.detectorPos2D, zeros(size(obj.detectorPos2D,1),1)]); + else + dm_det2d = distmatrix(obj.detectorPos2D, zeros); + end + dm_det3d = distmatrix(obj.detectorPos3D); + exp = round(log10(mean(dm_det3d(dm_det3d(:)>0)))) - round(log10(mean(dm_det2d(dm_det2d(:)>0)))); + scaleFactor = 10^exp; + obj.detectorPos2D = obj.detectorPos2D*scaleFactor; + end + + + + % ------------------------------------------------------- function Project_3D_to_2D(obj) + if all(obj.sourcePos2D(:)==0) + obj.sourcePos2D = []; + end + if all(obj.detectorPos2D(:)==0) + obj.detectorPos2D = []; + end if isempty(obj.sourcePos2D) && isempty(obj.detectorPos2D) if isempty(obj.landmarkPos3D) || ~obj.isValidLandmarkLabels() nSource = size(obj.sourcePos3D,1); @@ -220,6 +252,7 @@ function Project_3D_to_2D(obj) % obj.sourcePos2D = convert_optodepos_to_circlular_2D_pos(obj.sourcePos3D, T, norm_factor); obj.detectorPos2D = convert_optodepos_to_circlular_2D_pos(obj.detectorPos3D, T, norm_factor); + obj.ScaleProjection2D(); end end end diff --git a/Group/DataTree/AcquiredData/Snirf/SnirfClass.m b/Group/DataTree/AcquiredData/Snirf/SnirfClass.m index 679d730..b396242 100644 --- a/Group/DataTree/AcquiredData/Snirf/SnirfClass.m +++ b/Group/DataTree/AcquiredData/Snirf/SnirfClass.m @@ -1178,7 +1178,11 @@ function SetMetaDataTags(obj, val) p = length(obj2.GetTime()); q = length(obj.aux(ii).GetTime()); end - datamat(:,ii) = resample(obj.aux(ii).GetDataTimeSeries(), p, q); + m = resamplee(obj.aux(ii).GetDataTimeSeries(), p, q); + if isrow(m) + m = m'; + end + datamat = [datamat, m]; %#ok end end diff --git a/Group/DataTree/Version.txt b/Group/DataTree/Version.txt index 0e1f39b..250f359 100644 --- a/Group/DataTree/Version.txt +++ b/Group/DataTree/Version.txt @@ -1 +1 @@ -1.17.2 \ No newline at end of file +1.17.4 \ No newline at end of file diff --git a/ImgRecon/ImageRecon.m b/ImgRecon/ImageRecon.m index 0fe23bc..5bfde51 100644 --- a/ImgRecon/ImageRecon.m +++ b/ImgRecon/ImageRecon.m @@ -38,7 +38,7 @@ end % default time range -set(handles.time_range, 'String',num2str([5 10])); +set(handles.time_range, 'String',num2str([-5 10])); % default alpha (regularization) for brain only reconstruction set(handles.alpha_brainonly, 'String',1e-2); @@ -49,13 +49,17 @@ % default beta (regularization) for brain and scalp reconstruction set(handles.beta_brain_scalp, 'String',1e-2); -trange = abs(atlasViewer.dataTree.currElem.GetVar('trange')); -set(handles.time_range,'String',num2str(trange)); - - +% Set default time range +dcAvg = atlasViewer.dataTree.currElem.GetVar('dcAvg'); +trange = atlasViewer.dataTree.currElem.GetVar('trange'); +if ~isempty(dcAvg) + trange = [dcAvg.time(1), dcAvg.time(end)]; +end +set(handles.time_range,'String',sprintf('%0.3f %0.3f', ceiln(trange(1),-3), floorn(trange(2),-3))); err = 0; + % --------------------------------------------------------------------------------- function ImageRecon_OpeningFcn(hObject, ~, handles, varargin) global atlasViewer @@ -223,6 +227,7 @@ function image_recon_Callback(~, ~, handles) %%%% Error checking of subject data itself if ErrorCheck_Data(dcAvg, tHRF, cond, trange) < 0 + close(h); return end @@ -257,14 +262,15 @@ function image_recon_Callback(~, ~, handles) % put A matrix together and combine with extinction coefficients E = GetExtinctions([SD.Lambda(1) SD.Lambda(2)]); E = E/10; %convert from /cm to /mm E raws: wavelength, columns 1:HbO, 2:HbR - Amatrix = [squeeze(Adot(:,:,1))*E(1,1) squeeze(Adot(:,:,1))*E(1,2); - squeeze(Adot(:,:,2))*E(2,1) squeeze(Adot(:,:,2))*E(2,2)]; + Amatrix = [squeeze(Adot(:,:,1))*E(1,1), squeeze(Adot(:,:,1))*E(1,2); + squeeze(Adot(:,:,2))*E(2,1), squeeze(Adot(:,:,2))*E(2,2)]; alpha = str2num(get(handles.alpha_brainonly,'String')); [HbO, HbR, err] = hmrImageReconConc(yavgimg, [], alpha, Amatrix); if err==1 MenuBox('Error: Number of channels in measuremnt is not the same as in Adot.', 'Okay'); + close(h); return; end @@ -284,13 +290,11 @@ function image_recon_Callback(~, ~, handles) elseif value2 == 1 % brain and scalp reconstruction without short separation regression (Zhan2012) - % Adot = Adot(activeChIdxs,:,:); - % Adot = Adot(longSepChLst,:,:); - Adot = Adot(activeChLst_SDpairs,:,:); + Adot = Adot(activeChIdxs,:,:); + Adot = Adot(longSepChLst,:,:); - % Adot_scalp = Adot_scalp(activeChIdxs,:,:); - % Adot_scalp = Adot_scalp(longSepChLst,:,:); - Adot_scalp = Adot_scalp(activeChLst_SDpairs,:,:); + Adot_scalp = Adot_scalp(activeChIdxs,:,:); + Adot_scalp = Adot_scalp(longSepChLst,:,:); % get alpha and beta for regularization alpha = str2num(get(handles.alpha_brain_scalp,'String')); %#ok<*ST2NM> @@ -545,3 +549,5 @@ function plotHb_Callback(~, ~, handles) end err = 0; + + diff --git a/ImgRecon/hmrConc2OD.m b/ImgRecon/hmrConc2OD.m index 621ef28..5ec4047 100644 --- a/ImgRecon/hmrConc2OD.m +++ b/ImgRecon/hmrConc2OD.m @@ -20,34 +20,56 @@ % OUTPUTS: % dod: the change in OD (#time points x #channels) % - function dod = hmrConc2OD( dc, SD, ppf ) nWav = length(SD.Lambda); -ml = SD.MeasList; +ml0 = SD.MeasList; +nTpts = size(dc,1); + +dod0 = zeros(size(dc,1),size(ml0,1), size(dc,4)); if length(ppf)~=nWav errordlg('The length of PPF must match the number of wavelengths in SD.Lambda'); - dod = zeros(size(dod,1),size(ml,1)); + dod = dod0; return end -nTpts = size(dc,1); - e = GetExtinctions( SD.Lambda ); e = e(:,1:2) / 10; % convert from /cm to /mm -%einv = inv( e'*e )*e'; + +% dod0 assumes sorted channel order. So we have to reorder an unsorted ml +% to make sure its channel order matches the implicit order of dod0. +% After we generate dod0, we reorder it's columns to match the orgiginal +% channel order in ml0 +[ml, order] = sortMl(ml0); 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),:)); - try - dod(:,[idx1 idx2']) = (e * dc(:,1:2,idx)')' .* (ones(nTpts,1)*rho*ppf); - catch - d = 1; + dod0(:,[idx1, idx2']) = (e * dc(:,1:2,idx)')' .* (ones(nTpts,1)*rho*ppf); end -% dc(:,:,idx) = ( einv * (dod(:,[idx1 idx2'])./(ones(nTpts,1)*rho*ppf))' )'; +dod = dod0(:,order,:); + + + +% ---------------------------------------------------------- +function [ml, order] = sortMl(ml0) +order = zeros(size(ml0,1),1); + +k1 = find( ml0(:,4)==1 ); +k2 = find( ml0(:,4)==2 ); + +ml1 = ml0(k1,:); +ml2 = ml0(k2,:); + +ml1 = sortrows(ml1); +ml2 = sortrows(ml2); + +ml = [ml1; ml2]; + +for iM = 1:size(ml0,1) + order(iM) = find(ml(:,1)==ml0(iM,1) & ml(:,2)==ml0(iM,2) & ml(:,4)==ml0(iM,4)); end -%dc(:,3,:) = dc(:,1,:) + dc(:,2,:); + diff --git a/Probe/compatibleProbes.m b/Probe/compatibleProbes.m index 6683fde..f18b47b 100644 --- a/Probe/compatibleProbes.m +++ b/Probe/compatibleProbes.m @@ -31,7 +31,7 @@ if size(probe1.ml,1) ~= size(probe2.ml,1) return; end - if ~all(probe1.ml(:) == probe2.ml(:)) + if ~all(all(probe1.ml(:,[1,2,4]) == probe2.ml(:,[1,2,4]))) return; end end diff --git a/Probe/findUniqueProbeFiles.m b/Probe/findUniqueProbeFiles.m index e2773ad..e002908 100644 --- a/Probe/findUniqueProbeFiles.m +++ b/Probe/findUniqueProbeFiles.m @@ -66,11 +66,14 @@ probe2.ml = sort(probe2.ml); if ~isempty(probe1.ml) && ~isempty(probe2.ml) if size(probe1.ml,1) ~= size(probe2.ml,1) + b = false; return; end - if ~all(probe1.ml(:) == probe2.ml(:)) + if ~(all(all(probe1.ml(:,[1,2,4]) == probe2.ml(:,[1,2,4])))) + b = false; return; end + return else return; end diff --git a/Probe/initProbe.m b/Probe/initProbe.m index 56984ef..e645f89 100644 --- a/Probe/initProbe.m +++ b/Probe/initProbe.m @@ -2,7 +2,6 @@ probe = struct( ... 'name','probe', ... 'pathname',filesepStandard(pwd), ... - 'filename_to_save','', ... 'handles',struct( ... 'labels',[], ... 'circles',[], ... diff --git a/Probe/isPreRegisteredProbe.m b/Probe/isPreRegisteredProbe.m index 5f2d272..814d67e 100644 --- a/Probe/isPreRegisteredProbe.m +++ b/Probe/isPreRegisteredProbe.m @@ -79,12 +79,8 @@ return end if probeHasSpringRegistration(probe) - try refpts1.labels = probe.registration.al(:,2); refpts1.pos = probe.optpos(cell2array(probe.registration.al(:,1)), :); - catch - d = 1; - end else refpts1 = probe.registration.refpts; end diff --git a/Probe/pullProbeToHeadsurf.m b/Probe/pullProbeToHeadsurf.m index 735a275..0ee78fa 100644 --- a/Probe/pullProbeToHeadsurf.m +++ b/Probe/pullProbeToHeadsurf.m @@ -1,6 +1,6 @@ function probe = pullProbeToHeadsurf(probe, head) if isPreRegisteredProbe(probe)==2 - MenuBox('Probe already registered ot head surface'); + MenuBox('Probe already registered to head surface'); return end if isempty(probe.optpos_reg) diff --git a/Utils/Shared/Version.txt b/Utils/Shared/Version.txt index afa2b35..b9268da 100644 --- a/Utils/Shared/Version.txt +++ b/Utils/Shared/Version.txt @@ -1 +1 @@ -1.8.0 \ No newline at end of file +1.8.1 \ No newline at end of file diff --git a/Utils/Shared/ceiln.m b/Utils/Shared/ceiln.m new file mode 100644 index 0000000..e51755f --- /dev/null +++ b/Utils/Shared/ceiln.m @@ -0,0 +1,8 @@ +function val = ceiln(val0, n) +if ~exist('n','var') + n = 0; +end +pw = ceil(log10(abs(val0))); % Find order of magnitude of val. +res = 10^(pw-n-2); % Resolution to round to. +valtemp = val0*res; +val = ceil(valtemp)/res; % < change floor() to ceil(), for ceiling equivalent. diff --git a/Utils/Shared/floorn.m b/Utils/Shared/floorn.m new file mode 100644 index 0000000..a7f50ab --- /dev/null +++ b/Utils/Shared/floorn.m @@ -0,0 +1,9 @@ +function val = floorn(val0, n) +if ~exist('n','var') + n = 0; +end +pw = ceil(log10(abs(val0))); % Find order of magnitude of val. +res = 10^(pw-n-2); % Resolution to round to. +valtemp = val0*res; +val = floor(valtemp)/res; % < change floor() to ceil(), for ceiling equivalent. + diff --git a/Utils/Shared/getAppDir.m b/Utils/Shared/getAppDir.m index dbde5ab..eb7e337 100644 --- a/Utils/Shared/getAppDir.m +++ b/Utils/Shared/getAppDir.m @@ -75,7 +75,9 @@ else dirname = fileparts(which('Homer3.m')); end - +if ~ispathvalid(dirname) + dirname = filesepStandard(pwd); +end dirname(dirname=='\') = '/'; if dirname(end) ~= '/' dirname(end+1) = '/'; diff --git a/Utils/Shared/pretty_print_struct.m b/Utils/Shared/pretty_print_struct.m index bb87eb2..79fd0e4 100644 --- a/Utils/Shared/pretty_print_struct.m +++ b/Utils/Shared/pretty_print_struct.m @@ -1,4 +1,4 @@ -function pretty_print_struct(st, indent, option, logger) +function pretty_print_struct(st, indent, option) spaces = ''; if ~exist('st','var') || isempty(st) @@ -10,11 +10,6 @@ function pretty_print_struct(st, indent, option, logger) if ~exist('option','var') || isempty(option) option = 1; end -if ~exist('logger','var') - logger = []; -end -logger = InitLogger(logger); - if iswholenum(indent) spaces = blanks(indent); elseif ischar(indent) @@ -22,6 +17,7 @@ function pretty_print_struct(st, indent, option, logger) end if isstruct(st) || isobject(st) + fields = properties(st); s = evalc('disp(st)'); c = str2cell_fast(s, char(10)); for ii=1:length(c) diff --git a/Utils/Shared/resamplee.m b/Utils/Shared/resamplee.m new file mode 100644 index 0000000..6555e01 --- /dev/null +++ b/Utils/Shared/resamplee.m @@ -0,0 +1,10 @@ +function resig = resamplee(sig,upsample,downsample) +if upsample*downsample<2^31 + resig = resample(sig,upsample,downsample); +else + sig1half = sig(1:floor(length(sig)/2)); + sig2half = sig(floor(length(sig)/2):end); + resig1half = resamplee(sig1half, floor(upsample/2), length(sig1half)); + resig2half = resamplee(sig2half, upsample-floor(upsample/2), length(sig2half)); + resig = [resig1half(:); resig2half(:)]; +end diff --git a/Version.txt b/Version.txt index 928f623..7ad1be3 100644 --- a/Version.txt +++ b/Version.txt @@ -1 +1 @@ -2.47.0 \ No newline at end of file +2.48.1 \ No newline at end of file From 3a3ca9a6e731f0bd71756a92b4dd6019d33f3509 Mon Sep 17 00:00:00 2001 From: sreekanth kura Date: Mon, 20 Nov 2023 20:44:42 -0500 Subject: [PATCH 06/20] FIx a bug in probe design in AV --- Probe/isPreRegisteredProbe.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Probe/isPreRegisteredProbe.m b/Probe/isPreRegisteredProbe.m index 5f2d272..68e27ab 100644 --- a/Probe/isPreRegisteredProbe.m +++ b/Probe/isPreRegisteredProbe.m @@ -81,7 +81,7 @@ if probeHasSpringRegistration(probe) try refpts1.labels = probe.registration.al(:,2); - refpts1.pos = probe.optpos(cell2array(probe.registration.al(:,1)), :); + refpts1.pos = probe.optpos_reg(cell2array(probe.registration.al(:,1)), :); catch d = 1; end From 31da485409b17bce1b469ef7f899f1ee91de2318 Mon Sep 17 00:00:00 2001 From: Meryem Ayse Yucel <49535526+mayucel@users.noreply.github.com> Date: Thu, 30 Nov 2023 18:46:16 -0500 Subject: [PATCH 07/20] Update ImageRecon.m Brain and scalp option was referring two old variable names. This fix is replacing them with activeChLst_SDpairs variable. --- ImgRecon/ImageRecon.m | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/ImgRecon/ImageRecon.m b/ImgRecon/ImageRecon.m index 5bfde51..51f713f 100644 --- a/ImgRecon/ImageRecon.m +++ b/ImgRecon/ImageRecon.m @@ -290,11 +290,10 @@ function image_recon_Callback(~, ~, handles) elseif value2 == 1 % brain and scalp reconstruction without short separation regression (Zhan2012) - Adot = Adot(activeChIdxs,:,:); - Adot = Adot(longSepChLst,:,:); - Adot_scalp = Adot_scalp(activeChIdxs,:,:); - Adot_scalp = Adot_scalp(longSepChLst,:,:); + Adot = Adot(activeChLst_SDpairs,:,:); + Adot_scalp = Adot_scalp(activeChLst_SDpairs,:,:); + % get alpha and beta for regularization alpha = str2num(get(handles.alpha_brain_scalp,'String')); %#ok<*ST2NM> From 293e1aea9e16ac127eb58d52799a7abb75283e9f Mon Sep 17 00:00:00 2001 From: Meryem Ayse Yucel <49535526+mayucel@users.noreply.github.com> Date: Fri, 1 Dec 2023 09:02:32 -0500 Subject: [PATCH 08/20] Update showReducedMesh.m Previously, if the user opted for the original mesh, it was still reducing the mesh. Cleaned up to code to provide both options (reduced and orig mesh). --- Utils/showReducedMesh.m | 52 +++++++++++++++++++---------------------- 1 file changed, 24 insertions(+), 28 deletions(-) diff --git a/Utils/showReducedMesh.m b/Utils/showReducedMesh.m index 83b7344..fcd4bd2 100644 --- a/Utils/showReducedMesh.m +++ b/Utils/showReducedMesh.m @@ -8,39 +8,35 @@ fv = mesh_orig; fvn = mesh_reduced; -mesh = fvn; + + +hf = figure; +% Plot the original +subplot(1,2,1) +h=trisurf( fv.faces, fv.vertices(:,1), fv.vertices(:,2), fv.vertices(:,3) ); +title( sprintf(' Original with %d faces', size(fv.faces,1)) ) +set(h,'linestyle','none') +light + +% Plot the reduced +subplot(1,2,2) +h=trisurf( fvn.faces, fvn.vertices(:,1), fvn.vertices(:,2), fvn.vertices(:,3) ); +title( sprintf(' New with %d faces', size(fvn.faces,1)) ) +set(h,'linestyle','none') +light q = MenuBox( sprintf('The original pial surface has %d faces. It is recommended that this be less than 40,000 for best performance.\nShall I reduce this?',size(fv.faces,1)),{'Yes','No'}); + if q==1 - hf = figure; - - % Plot the original - subplot(1,2,1) - h=trisurf( fv.faces, fv.vertices(:,1), fv.vertices(:,2), fv.vertices(:,3) ); - title( sprintf(' Original with %d faces', size(fv.faces,1)) ) - set(h,'linestyle','none') - light - - % Plot the reduced - subplot(1,2,2) - h=trisurf( fvn.faces, fvn.vertices(:,1), fvn.vertices(:,2), fvn.vertices(:,3) ); - title( sprintf(' New with %d faces', size(fvn.faces,1)) ) - set(h,'linestyle','none') - light - - q = MenuBox('Accept this?',{'Yes','No'}); - if q==2 - q = MenuBox('Proceed with the original mesh? This could take a long time.',{'Yes','No'}); - if q==1 - mesh = fv; - end - end - - close(hf); + mesh = fvn; +else + mesh = fv; end +close(hf); + + nodeX = mesh.vertices; -elem = mesh.faces; +elem = mesh.faces; nNode = size(nodeX,1); - From 3e49e7ae9d3755e1fb13ff5fe02708aefb8dc728 Mon Sep 17 00:00:00 2001 From: sreekanth kura Date: Mon, 4 Dec 2023 11:03:42 -0500 Subject: [PATCH 09/20] Fix isProbeAlreadyRegistered issue --- AtlasViewerGUI.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/AtlasViewerGUI.m b/AtlasViewerGUI.m index d796a04..e6f3f37 100644 --- a/AtlasViewerGUI.m +++ b/AtlasViewerGUI.m @@ -653,7 +653,7 @@ function pushbuttonRegisterProbeToSurface_Callback(hObject, eventdata, handles) refpts = set_eeg_active_pts(refpts, [], false); % Finish registration -if isPreRegisteredProbe(probe, refpts) +if 0 %isPreRegisteredProbe(probe, refpts) % check if probe is already on the surface probe_on_surface = true; From cb66adaba2253f18f07c5bf09c5ba765f9cb1727 Mon Sep 17 00:00:00 2001 From: sreekanth kura Date: Tue, 30 Jan 2024 14:01:11 -0500 Subject: [PATCH 10/20] update code to add transformation matrix from mcspace to mni space after importing a new MRI anatomy --- Labelssurf/getLabelssurf.m | 14 ++++++++++---- Labelssurf/saveLabelSurf2Vol.m | 20 ++++++++++++++++++++ Utils/fs2viewer/ImportMriAnatomy.m | 1 + 3 files changed, 31 insertions(+), 4 deletions(-) create mode 100644 Labelssurf/saveLabelSurf2Vol.m diff --git a/Labelssurf/getLabelssurf.m b/Labelssurf/getLabelssurf.m index 578eca9..98606b7 100644 --- a/Labelssurf/getLabelssurf.m +++ b/Labelssurf/getLabelssurf.m @@ -19,6 +19,12 @@ end dirname = [dirname0 'anatomical/']; +if exist([dirname 'labelssurf2vol.txt'],'file') + T_2vol = load([dirname 'labelssurf2vol.txt'],'ascii'); + labelssurf.T_2vol = T_2vol; +end + + if exist([dirname 'labelssurf.mat'],'file') load([dirname 'labelssurf.mat'],'-mat'); else @@ -33,9 +39,9 @@ end load([dirname 'labelssurf.mat']); -if exist([dirname 'labelssurf2vol.txt'],'file') - T_2vol = load([dirname 'labelssurf2vol.txt'],'ascii'); -end +% if exist([dirname 'labelssurf2vol.txt'],'file') +% T_2vol = load([dirname 'labelssurf2vol.txt'],'ascii'); +% end fv = struct('vertices',[],'faces',[]); idxL = []; @@ -65,7 +71,7 @@ labelssurf.colormaps = cm; labelssurf.colormapsIdx = 4; labelssurf.idxL = idxL; -labelssurf.T_2vol = T_2vol; +% labelssurf.T_2vol = T_2vol; if ~labelssurf.isempty(labelssurf) diff --git a/Labelssurf/saveLabelSurf2Vol.m b/Labelssurf/saveLabelSurf2Vol.m new file mode 100644 index 0000000..0b2bbd4 --- /dev/null +++ b/Labelssurf/saveLabelSurf2Vol.m @@ -0,0 +1,20 @@ +function saveLabelSurf2Vol(headvol) +% Author:Sreekanth Kura (skura@bu.edu) + +% I am assuming imported anatomy will be saved in RAS orientation. +S = size(headvol.img); +T = headvol.T_2ras; +Tras_2vol = T; +for u = 1:3 + if sum(Tras_2vol(u,1:3)) == -1 + Tras_2vol(u,4) = S(u); + else + Tras_2vol(u,4) = 0; + end +end + +Tras_mni = (Tras_2vol'*T')'; +Tmni_ras = inv(Tras_mni); +dirname = [headvol.pathname, '/anatomical/']; +filename = [dirname 'labelssurf2vol.txt']; +writematrix(Tmni_ras,filename,'Delimiter',' ') \ No newline at end of file diff --git a/Utils/fs2viewer/ImportMriAnatomy.m b/Utils/fs2viewer/ImportMriAnatomy.m index a0d2d7a..70c18a3 100644 --- a/Utils/fs2viewer/ImportMriAnatomy.m +++ b/Utils/fs2viewer/ImportMriAnatomy.m @@ -356,6 +356,7 @@ function ImportMriAnatomy_CloseRequestFcn(~, ~, handles) saveHeadvol(headvol); saveHeadsurf(headsurf); savePialsurf(pialsurf); + saveLabelSurf2Vol(headvol); catch ME From ff0cfac11a8a285aeec7aa8d590502959804f45e Mon Sep 17 00:00:00 2001 From: sreekanth kura Date: Thu, 8 Feb 2024 15:48:04 -0500 Subject: [PATCH 11/20] Fix merging error --- Probe/isPreRegisteredProbe.m | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Probe/isPreRegisteredProbe.m b/Probe/isPreRegisteredProbe.m index 9389456..bf5742c 100644 --- a/Probe/isPreRegisteredProbe.m +++ b/Probe/isPreRegisteredProbe.m @@ -81,9 +81,9 @@ if probeHasSpringRegistration(probe) refpts1.labels = probe.registration.al(:,2); refpts1.pos = probe.optpos_reg(cell2array(probe.registration.al(:,1)), :); - catch - d = 1; - end + % catch + % d = 1; + % end else refpts1 = probe.registration.refpts; end From 86324b0babe683869747593902a914c8652de308 Mon Sep 17 00:00:00 2001 From: jayd1860 Date: Sun, 18 Feb 2024 23:06:00 -0500 Subject: [PATCH 12/20] -- Add unit test script to test DataFilesClass with simulated dataset --- Group/UnitTests/testDataFilesClass.m | 109 +++++++++++++++++++++++++++ 1 file changed, 109 insertions(+) create mode 100644 Group/UnitTests/testDataFilesClass.m diff --git a/Group/UnitTests/testDataFilesClass.m b/Group/UnitTests/testDataFilesClass.m new file mode 100644 index 0000000..b23dc4d --- /dev/null +++ b/Group/UnitTests/testDataFilesClass.m @@ -0,0 +1,109 @@ +function [files0, files1, files2] = testDataFilesClass(options) +% Syntax: +% [files0, files1, files2] = testDataFilesClass() +% [files0, files1, files2] = testDataFilesClass(options) +% +% Examples: +% [files0, files1, files2] = testDataFilesClass() +% [files0, files1, files2] = testDataFilesClass('create') +% [files0, files1, files2] = testDataFilesClass('create:remove') +% +% +if ~exist('options','var') + options = ''; +end +datasetDir = './group1'; + +if includes(options,'create') + if ~exist(datasetDir,'dir') + fprintf('mkdir(''%s'');\n', datasetDir) + mkdir(datasetDir); + end + fprintf('cd(''%s'');\n', datasetDir) + cd(datasetDir) +end +dirs = dir('sub-*'); +if isempty(dirs) + fprintf('generateSimDataset(pwd, 5, 3, 4);\n') + generateSimDataset(pwd, 5, 3, 4); +end +fprintf('\n') +pause(2) + + +if exist('./derivatives','dir') + fprintf('rmdir(''./derivatives'',''s'')\n'); + rmdir('./derivatives','s') +end +fprintf('tic; files0 = DataFilesClass(pwd, ''snirf''); toc;\n'); +tic; files0 = DataFilesClass(pwd, 'snirf'); toc %#ok<*NASGU> +fprintf('\n') +pause(2) + +% Update files +filesUpdated = { + './sub-02/ses-01/sub-02_ses-01_run02.snirf'; + './sub-02/ses-03/sub-02_ses-03_run01.snirf'; + './sub-03/ses-02/sub-03_ses-02_run02.snirf'; + './sub-03/ses-02/sub-03_ses-02_run04.snirf'; + './sub-05/ses-01/sub-05_ses-01_run04.snirf'; + }; +s = repmat(SnirfClass(),1,length(filesUpdated)); +for ii = 1:length(filesUpdated) + fprintf('Update file %s\n', filesUpdated{ii}) + s(ii) = SnirfClass(filesUpdated{ii}); + s(ii).Save(); +end +fprintf('\n') +pause(2) + + +% Added files +filesAdded = { + './sub-01/ses-02/sub-01_ses-02_run10.snirf'; + './sub-03/ses-02/sub-04_ses-02_run14.snirf'; + './sub-04/ses-02/sub-04_ses-02_run11.snirf'; + './sub-04/ses-03/sub-04_ses-03_run12.snirf'; + }; +for ii = 1:length(filesAdded) + fprintf('Added file %s\n', filesAdded{ii}) + s = SnirfClass('./sub-01/ses-01/sub-01_ses-01_run01.snirf'); + s.Save(filesAdded{ii}); +end +fprintf('\n') +pause(2) + + +% Deleted files +filesDeleted = { + './sub-02/ses-03/sub-02_ses-03_run02.snirf'; + './sub-04/ses-02/sub-04_ses-02_run01.snirf'; + './sub-05/ses-02/sub-05_ses-02_run03.snirf'; + }; +for ii = 1:length(filesDeleted) + fprintf('Delete file %s\n', filesDeleted{ii}) + delete(filesDeleted{ii}) +end +fprintf('\n') +pause(2) + + +fprintf('tic; files1 = DataFilesClass(pwd, ''snirf''); toc;\n'); +tic; files1 = DataFilesClass(pwd, 'snirf'); toc +fprintf('\n') +pause(2) + +fprintf('tic; files2 = DataFilesClass(pwd, ''snirf''); toc;\n'); +tic; files2 = DataFilesClass(pwd, 'snirf'); toc +fprintf('\n') +pause(2) + +if includes(options,'remove') + fprintf('cd(''..'');\n') + cd('..') + fprintf('rmdir(''%s'',''s'');\n', datasetDir) + rmdir(datasetDir,'s') + fprintf('\n') +end + + From 85b38413d9a66b5c2de7bbb6fc7afb4c13728ad4 Mon Sep 17 00:00:00 2001 From: jayd1860 Date: Sun, 18 Feb 2024 23:13:40 -0500 Subject: [PATCH 13/20] v2.49.0 * DataTree, v1.19.0 -- Building on last commit to improve DataTree launch performance by caching DataFilesClass object in ./derivatives/homer/DatasetFiles_snirf.mat - this time add smart loading of ./derivatives/homer/DatasetFiles_snirf.mat containing cached DataFilesClass object which takes a long time to build. Rather than just crudely rebuilding DataFilesClass object from scratch if anything in the current files differs from saved data; rebuild ONLY those parts of the dataset that have changes and add those changes to saved object. * DataTree, v1.18.1 -- Changes to improve startup performance when launching Homer. It takes a lot of time to build and check the set of valid data files out of which the dataTree is built. changes to BUNPC/development. It will save the object with the vetted dataset files list (that's what slows things down) in the BIDS derived folder derivatives / homer / DatasetFiles.mat It will do it once if the saved file isn't there OR if something about the dataset has changed since last save - i.e. data file added, deleted, moved or edited. Otherwise it quickly loads and copies it which is a lot faster. -- Fix issue with incorrect folder structure report in DataFilesClass.m and incorrect counting of number of files and folders. Test case is TaskA -- Fix ProbeClass.IsValid validation not checking mandatory dimensions of sourcePos and detectorPos os having minimum of 2 and 3 columns. * DataTree, 1.18.0 -- Fix for "Stim Rejection not working for auto/manual time exclusion" --- .../AcquiredData/DataFiles/DataFilesClass.m | 323 +++++++++++++++--- .../AcquiredData/DataFiles/FileClass.m | 70 +++- .../DataTree/AcquiredData/Snirf/ProbeClass.m | 20 ++ .../DataTree/AcquiredData/Snirf/SnirfClass.m | 20 ++ Group/DataTree/AcquiredData/Snirf/StimClass.m | 26 +- Group/DataTree/Version.txt | 2 +- Version.txt | 2 +- 7 files changed, 410 insertions(+), 53 deletions(-) diff --git a/Group/DataTree/AcquiredData/DataFiles/DataFilesClass.m b/Group/DataTree/AcquiredData/DataFiles/DataFilesClass.m index a118540..1d08907 100644 --- a/Group/DataTree/AcquiredData/DataFiles/DataFilesClass.m +++ b/Group/DataTree/AcquiredData/DataFiles/DataFilesClass.m @@ -16,6 +16,8 @@ properties (Access = private) lookupTable excludedFolders + savedFilename + changed end methods @@ -30,6 +32,8 @@ obj.err = -1; obj.errmsg = {}; obj.dirFormats = struct('type',0, 'choices',{{}}); + obj.lookupTable = []; + obj.changed = false; global logger global cfg @@ -38,7 +42,6 @@ cfg = InitConfig(cfg); obj.logger = logger; - obj.lookupTable = []; obj.excludedFolders = {}; skipconfigfile = false; @@ -76,7 +79,7 @@ obj.rootdir = filesepStandard(obj.rootdir,'full'); % Configuration parameters - obj.config = struct('RegressionTestActive','', 'AskToFixNameConflicts',1, 'DerivedDataDir',''); + obj.config = struct('RegressionTestActive','', 'AskToFixNameConflicts',1, 'DerivedDataDir','', 'DerivedDataRootDir',''); if skipconfigfile==false str = cfg.GetValue('Regression Test Active'); if strcmp(str,'true') @@ -110,108 +113,287 @@ [obj.rootdir, 'imagerecon']; [obj.rootdir, 'anatomical']; }; + obj.config.DerivedDataRootDir = [obj.rootdir, obj.config.DerivedDataDir, '/', f, '/']; + obj.savedFilename = [obj.config.DerivedDataRootDir, 'DatasetFiles', '_', obj.filetype, '.mat']; if nargin==0 return; end obj.err = 0; - - obj.InitFolderFromats(); + obj.InitFolderFormats(); obj.GetDataSet(); + obj.SaveDataSet(); end % ----------------------------------------------------------------------------------- - function GetDataSet(obj) + function Compare(obj) + if obj.IsEmpty() + return + end + obj2 = DataFilesClass(); + obj2.rootdir = obj.rootdir; + obj2.filetype = obj.filetype; + obj2.dirFormats = obj.dirFormats; + obj2.GetDataSet('skipsaveddata') + + %%%% Create list of all the saved names for sorting later + filesSaved = cell(length(obj.files),1); + filesSavedErr = cell(length(obj.filesErr),1); + for ii = 1:length(obj.files) + filesSaved{ii} = obj.files(ii).name; + end + for ii = 1:length(obj.filesErr) + filesSavedErr{ii} = obj.filesErr(ii).name; + end + + %%%% Move any files from saved error list that have been updated to nonerror list + N = length(obj.filesErr); + for ii = N:-1:1 + for jj = 1:length(obj2.files) + if strcmp(obj.filesErr(ii).name, obj2.files(jj).name) + if obj.filesErr(ii).datenum ~= obj2.files(jj).datenum + obj.filesErr(ii) = []; + filesSavedErr{ii} = []; + obj.files(end+1) = obj2.files(jj); + if ~obj2.files(jj).isdir + obj.nfiles = obj.nfiles+1; + end + filesSaved{end+1} = obj2.files(jj).name; %#ok + obj.logger.Write('%d. Mismatch in dates for saved data in %s. Will re-validate this file.\n', jj, obj.filesErr(ii).name); + obj.changed = true; + end + break + end + end + end + + + %%%% Now check saved object nonerror list + N = length(obj.files); + for ii = N:-1:1 + found = false; + for jj = 1:length(obj2.files) + if strcmp(obj.files(ii).name, obj2.files(jj).name) + found = true; + if obj.files(ii).datenum ~= obj2.files(jj).datenum + obj.files(ii) = obj2.files(jj).copy(); + obj.logger.Write('%d. Mismatch in dates for saved data in %s. Will re-validate this file.\n', jj, obj.files(ii).name); + obj.changed = true; + end + break + end + end + if ~found + obj.logger.Write('%d. Deleting file %s from saved data.\n', ii, obj.files(ii).name); + if ~obj.files(ii).isdir + obj.nfiles = obj.nfiles-1; + end + obj.files(ii) = []; + filesSaved(ii) = []; + obj.changed = true; + end + end + + % We don't need to check if files were added if there's no + % change and the number of current files found is same as saved + % number of files. Exist early. + if obj.changed == false + if (length(obj.files) + length(obj.filesErr)) == length(obj2.files) + return + end + end + + + %%%% Find which files were added to current dataset and add + %%%% them to saved object + for ii = 1:length(obj2.files) + found = false; + for jj = 1:length(obj.files) + if strcmp(obj2.files(ii).name, obj.files(jj).name) + found = true; + break + end + end + + % If the current file filesCurr(ii) is not found in saved dataset then it's a new file. + % Add new file from current dataset to saved dataset + if ~found + obj.files(end+1) = obj2.files(ii); + if ~obj2.files(ii).isdir + obj.nfiles = obj.nfiles+1; + end + filesSaved{end+1} = obj2.files(ii).name; %#ok + obj.logger.Write('%d. Adding new file %s to saved data.\n', jj, obj2.files(ii).name); + obj.changed = true; + end + end + + + % Sort non-error files + if obj.changed + [~, order] = sort(filesSaved); + obj.files = obj.files(order); + obj.logger.Write('\n'); + end + + end + + + + % ----------------------------------------------------------------------------------- + function SaveDataSet(obj) + if ~obj.changed + return + end + obj.changed = false; + obj.logger.Write('Saving DataFilesClass object in %s\n', obj.savedFilename); + if ~exist(obj.config.DerivedDataRootDir, 'dir') + mkdir(obj.config.DerivedDataRootDir) + end + save(obj.savedFilename, '-mat', 'obj'); + end + + + + % ----------------------------------------------------------------------------------- + function b = LoadDataSet(obj) + b = false; + if ~exist(obj.savedFilename, 'file') + return + end + dataset = load(obj.savedFilename); + dataset.obj.Compare(); + if isempty(dataset.obj.files) && isempty(dataset.obj.filesErr) + return + end + obj.logger.Write('Loaded saved DataFilesClass object: %s\n', obj.savedFilename); + obj.Copy(dataset.obj); + + obj.ErrorCheck() + obj.ErrorCheckName() + b = true; + end + + + + % ----------------------------------------------------------------------------------- + function GetDataSet(obj, skipsaveddata) + if ~exist('skipsaveddata','var') + skipsaveddata = ''; + end + if ~strcmp(skipsaveddata, 'skipsaveddata') + if obj.LoadDataSet() + return + end + end + + obj.changed = true; + + mode = ''; + if obj.dirFormats.type == 0 + iFormats = 1:length(obj.dirFormats.choices); + else + iFormats = obj.dirFormats.type; + mode = 'noerrorcheck'; + end if exist(obj.rootdir, 'dir')~=7 error('Invalid subject folder: ''%s''', obj.rootdir); end - - for ii = 1:length(obj.dirFormats.choices) + for ii = iFormats obj.InitLookupTable(); obj.FindDataSet(ii); + if strcmp(mode,'noerrorcheck') + continue + end + if isempty(obj.files) + continue + end + % Remove any files that cannot pass the basic test of loading % its data obj.ErrorCheck(); if ~isempty(obj.files) break end + end + if strcmp(mode,'noerrorcheck') + return end - obj.ErrorCheckFinal(); end % ----------------------------------------------------------------------------------- - function InitFolderFromats(obj) + function InitFolderFormats(obj) obj.dirFormats.choices = { %%%% 1. Flat #1 { - ['*_run*.', obj.filetype]; + ['*_run*.', obj.filetype]; } %%%% 2. Flat #2 { - ['*.', obj.filetype]; + ['*.', obj.filetype]; } %%%% 3. Subjects { - '*'; - ['*.', obj.filetype]; + '*'; + ['*.', obj.filetype]; } %%%% 4. BIDS #1, sub-