Update setpaths for generalizability

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
@@ -651,29 +653,12 @@ 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;
-    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 +1266,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 +1495,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 +1575,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 +1648,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 +1668,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 +1940,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 +2227,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 +2245,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 +2277,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 +2355,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 +2365,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 +2456,7 @@ function popupmenuImageDisplay_Callback(hObject, ~, handles)
 
 
 % --------------------------------------------------------------------
-function popupmenuImageDisplay_CreateFcn(hObject, eventdata, handles)
+function popupmenuImageDisplay_CreateFcn(hObject, ~, ~)
 global popupmenuorder;
 
 popupmenuorder = struct(...

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)

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