registration optimizations

- we now can have sessions with different mouse_names/dates!
- minor memory management improvements (getBatchSize())
- speeding up movie registration step on GPU (really using
batchProcessing now, was not previously)
- significantly speeding up loading the tiff files from disk (local SSD)
using memmapfile()

build_ops3.m

@@ -1,25 +1,54 @@
 function ops = build_ops3(db, ops)
 
 % ops = db;
-ops = addfields(ops, db);
-
-for k = 1:length(db.expts)
-    ops.SubDirs{k}    = num2str(db.expts(k));
-end
-
-ops.RootDir = fullfile(ops.RootStorage, ops.mouse_name, ops.date);
-
-% build file list
-ops.fsroot = [];
-for j = 1:length(ops.SubDirs)
-    ops.fsroot{j} = dir(fullfile(ops.RootDir, ops.SubDirs{j}, '*.tif'));
-    for k = 1:length(ops.fsroot{j})
-        ops.fsroot{j}(k).name = fullfile(ops.RootDir, ops.SubDirs{j}, ops.fsroot{j}(k).name);
+if ~iscell(db.mouse_name) 
+    % this is the usual case where we have a simple single session recording
+    ops = addfields(ops, db);
+
+    for k = 1:length(db.expts)
+        ops.SubDirs{k}    = num2str(db.expts(k));
+    end
+
+    ops.RootDir = fullfile(ops.RootStorage, ops.mouse_name, ops.date);
+
+    % build file list
+    ops.fsroot = [];
+    for j = 1:length(ops.SubDirs)
+        ops.fsroot{j} = dir(fullfile(ops.RootDir, ops.SubDirs{j}, '*.tif'));
+        for k = 1:length(ops.fsroot{j})
+            ops.fsroot{j}(k).name = fullfile(ops.RootDir, ops.SubDirs{j}, ops.fsroot{j}(k).name);
+        end
+    end
+else
+    % here we might have multiple sessions, which we want to be analyzed
+    % together (exactly the same FOV)
+    nSessions = length(db.mouse_name);
+    % a backwards compatible version of db
+    dbCompat = db;
+    dbCompat.mouse_name = db.mouse_name{1};
+    dbCompat.date = db.date{1};
+    dbCompat.expts = cell2mat(db.expts(:)');
+    ops = addfields(ops, dbCompat);
+    ops.db_orig = db;
+
+    ops.fsroot = cell(0);
+    ops.SubDirs = cell(0);
+    for iSession = 1:nSessions
+        ops.RootDir = fullfile(ops.RootStorage, db.mouse_name{iSession}, db.date{iSession});
+        for iExp = 1:length(db.expts{iSession})
+            ops.SubDirs{end+1} = num2str(db.expts{iSession}(iExp));
+            ops.fsroot{end+1} = dir(fullfile(ops.RootDir, ops.SubDirs{end}, '*.tif'));
+            for iFile = 1:length(ops.fsroot{end})
+                ops.fsroot{end}(iFile).name = fullfile(ops.RootDir, ops.SubDirs{end}, ops.fsroot{end}(iFile).name);
+            end
+        end
     end
+    % this line to be backward compatible (just in case)
+    ops.RootDir = fullfile(ops.RootStorage, db.mouse_name{1}, db.date{1});
 end
 
-try    
-     % MK code for automatically determining number of planes and channels
+try
+    % MK code for automatically determining number of planes and channels
     [~, header] = loadFramesBuff(ops.fsroot{1}(1).name, 1, 1, 1);
 
     hh=header{1};
@@ -50,7 +79,7 @@
     ops.planesToProcess = 1:ops.nplanes;
 else
     % planesToProcess is not working right now
-    ops.planesToProcess = 1:ops.nplanes; 
+    ops.planesToProcess = 1:ops.nplanes;
 end
 
 CharSubDirs = '';
@@ -60,4 +89,4 @@
 CharSubDirs = CharSubDirs(1:end-1);
 
 ops.ResultsSavePath = sprintf('%s//%s//%s//%s//', ops.ResultsSavePath, ops.mouse_name, ops.date, ...
-        CharSubDirs);
+    CharSubDirs);

configFiles/make_db_example.m

@@ -15,10 +15,18 @@
 db(i).date          = '2015-12-19';
 db(i).expts         = [4];
 
+% example of datasets, which consist of several sessions - use cell arrays
+% will be treated as subsets of experiment with the same FOV, with
+% different names/dates (for one reason or another), analyzed together
+i = i+1;
+db(i).mouse_name    = {'MK020', 'M150416_MK020'};
+db(i).date          = {'2015-07-30', '2015-07-30'};
+db(i).expts         = {[2010 2107], [1 2 3]};
+
 % example extra entries
 % db(i).AlignToRedChannel= 1;
 % db(i).BiDiPhase        = 0; % adjust the relative phase of consecutive lines
 % db(i).nSVD             = 1000; % will overwrite the default, only for this dataset
 % db(i).comments      = 'this was an adaptation experiment';
 % db(i).expred        = [4]; % say one block which had a red channel 
-% db(i).nchannels_red = 2; % how many channels did the red block have in total (assumes red is last)
+% db(i).nchannels_red = 2; % how many channels did the red block have in total (assumes red is last)

registration/getBatchSize.m

@@ -0,0 +1,8 @@
+function batchSize = getBatchSize(nPixels)
+
+g = gpuDevice;
+batchSize = 2^(floor(log2(g.AvailableMemory))-6)/2^ceil(log2(nPixels));
+
+% The calculation was deducted from the following examples
+% batchSize = 2^25/2^ceil(log2(nPixels)); % works well on GTX 970 (which has 4 GB memory)
+% batchSize = 2^23/2^ceil(log2(nPixels)); % works well on GTX 560 Ti (which has 1 GB memory)

registration/reg2P.m

@@ -145,15 +145,16 @@
     for j = 1:length(fs{k})
         iplane0 = mod(iplane0-1, numPlanes) + 1;
 
-        nFr = nFramesTiff(fs{k}(j).name);
         if red_align
             %                 ichanset = [nchannels*(iplane0-1)+rchannel; nFr; nchannels];
-            ichanset = [rchannel; nFr; nchannels];
+            ichanset = [rchannel; NaN; nchannels];
         else
-            ichanset = [ichannel; nFr; nchannels];
+            % using NaN instead of nFrames, we want to ask for nFrames
+            % after copying the file locally (way faster)
+            ichanset = [ichannel; NaN; nchannels];
         end
         data = loadFramesBuff(fs{k}(j).name, ichanset(1), ichanset(2), ichanset(3), ops.temp_tiff);
-
+        nFr = size(data, 3)*nchannels;
 
         if BiDiPhase
             yrange = 2:2:Ly;
@@ -190,17 +191,20 @@
             % if aligning by the red channel, data needs to be reloaded as the
             % green channel
             if red_align
+                nFr = nFramesTiff(fs{k}(j).name);
                 if mod(nFr, nchannels) ~= 0
                     fprintf('  WARNING: number of frames in tiff (%d) is NOT a multiple of number of channels!\n', j);
                 end
                 ichanset = [ichannel; nFr; nchannels];
-%                 data = img.loadFrames(fs{k}(j).name, ichanset(1), ichanset(2), ichanset(3));                
-                data = loadFramesBuff(ops.temp_tiff, ichanset(1), ichanset(2), ichanset(3));               
+%                 data = img.loadFrames(fs{k}(j).name, ichanset(1), ichanset(2), ichanset(3));           
+                % providing the last argument ops.temp_tiff to let the
+                % function know that this is already a local file (will load faster)
+                data = loadFramesBuff(ops.temp_tiff, ichanset(1), ichanset(2), ichanset(3), ops.temp_tiff);               
             end
 
             ix0 = 0;
             Nbatch = 1000;
-            dreg = zeros(size(data), class(data));            
+            dreg = zeros(size(data), 'like', data);            
             while ix0<size(data,3)
                 indxr = ix0 + (1:Nbatch);
                 indxr(indxr>size(data,3)) = [];

registration/register_movie.m

@@ -2,36 +2,58 @@
 
 orig_class = class(data);
 
-if ops.useGPU
-    data = gpuArray(single(data));
-end
-[Ly, Lx, NT] = size(data);
-
-Ny = ifftshift([-fix(Ly/2):ceil(Ly/2)-1]);
-Nx = ifftshift([-fix(Lx/2):ceil(Lx/2)-1]);
-[Nx,Ny] = meshgrid(Nx,Ny);
-Nx = Nx / Lx;
-Ny = Ny / Ly;
-
-if ops.useGPU
-    dreg = gpuArray.zeros(size(data), orig_class);
-else
-    dreg = zeros(size(data), orig_class);
-end
-
-
-if ops.useGPU
-    ds = gpuArray(ds);
-    Nx = gpuArray(single(Nx));
-    Ny = gpuArray(single(Ny));
-end
-
-for i = 1:NT
-    dph         = 2*pi*(ds(i,1)*Ny + ds(i,2)*Nx);
-    fdata       = fft2(single(data(:,:,i)));
-    dreg(:,:,i) = real(ifft2(fdata .* exp(1i * dph)));
-end
-
-if ops.useGPU
-    dreg = gather_try(dreg);
+[h, w, nFrames] = size(data);
+nFramesPerBatch = getBatchSize(w*h);
+nBatches = ceil(nFrames/nFramesPerBatch);
+startFrame = 1:nFramesPerBatch:nFrames;
+endFrame = min(startFrame+nFramesPerBatch-1, nFrames);
+dreg = zeros(size(data), orig_class);
+
+
+for iBatch = 1:nBatches
+    idx = startFrame(iBatch):endFrame(iBatch);
+    if ops.useGPU
+        dataBatch = gpuArray(single(data(:,:,idx)));
+    else
+        dataBatch = data(:,:,idx);
+    end
+    [Ly, Lx, NT] = size(dataBatch);
+
+    Ny = ifftshift([-fix(Ly/2):ceil(Ly/2)-1]);
+    Nx = ifftshift([-fix(Lx/2):ceil(Lx/2)-1]);
+    [Nx,Ny] = meshgrid(Nx,Ny);
+    Nx = Nx / Lx;
+    Ny = Ny / Ly;
+
+    if ops.useGPU
+        dregBatch = gpuArray.zeros(size(dataBatch), orig_class);
+    else
+        dregBatch = zeros(size(dataBatch), orig_class);
+    end
+
+    if ops.useGPU
+        dsBatch = gpuArray(permute(ds(idx, :), [3, 2, 1]));
+        Nx = gpuArray(single(Nx));
+        Ny = gpuArray(single(Ny));
+    else
+        dsBatch = ds(idx, :);
+    end
+
+    if ops.useGPU % do it batch-by-batch
+        dph         = 2*pi*(bsxfun(@times, dsBatch(1,1,:), Ny) + ...
+            bsxfun(@times, dsBatch(:,2,:), Nx));
+        fdata       = fft2(dataBatch);
+        dregBatch = real(ifft2(fdata .* exp(1i * dph)));
+    else % do it frame-by-frame
+        for i = 1:NT
+            dph         = 2*pi*(dsBatch(i,1)*Ny + dsBatch(i,2)*Nx);
+            fdata       = fft2(single(dataBatch(:,:,i)));
+            dregBatch(:,:,i) = real(ifft2(fdata .* exp(1i * dph)));
+        end
+    end
+
+    if ops.useGPU
+        dregBatch = gather_try(dregBatch);
+    end
+    dreg(:,:,idx) = dregBatch;
 end

registration/registration_offsets.m

@@ -55,7 +55,7 @@
   cfRefImg = cfRefImg./(eps0 + abs(cfRefImg)).*fhg;
 end
 if useGPU
-  batchSize = 2^25/2^ceil(log2(lyus*lxus)); % works well on GTX 970
+  batchSize = getBatchSize(lyus*lxus);
   maskMul = gpuArray(maskMul);
   maskOffset = gpuArray(maskOffset);
   cfRefImg = gpuArray(cfRefImg);