Fix correctness issues in readers, dF/F, deconvolution and FluFinder detection

.github/workflows/test-code.yml

@@ -21,6 +21,7 @@ concurrency:
 permissions:
   contents: write
   checks: write
+  pull-requests: write
   security-events: write
 
 jobs:

code/+nansen/+module/+twophoton/+autosegmentation/+flufinder/+detect/findUniqueRoisFromComponents.m

@@ -58,17 +58,28 @@
     % Boolean for all remaining components that wasn't taken care of yet
     remaining = true(numel(S), 1);
 
-    % Create a "sum projection" image of all components.
-    uniquePixelList = unique(allPixelIndices);
-
+    % Create a "sum projection" image of all components: the number of
+    % components each pixel belongs to. accumarray counts every pixel,
+    % including the last unique one, and copes with zero or one pixel
+    % (histcounts needs at least two bin edges).
     componentImage = zeros(imageSize);
-    % Using a histogram here is faster than a for loop.      % (Maybe only if there are 100ks of components...)
-    [N,E] = histcounts(allPixelIndices, uniquePixelList);
-    componentImage(E(1:end-1)) = N;
+    if ~isempty(allPixelIndices)
+        componentImage(:) = accumarray(allPixelIndices(:), 1, [numel(componentImage), 1]);
+    end
 
     summary.ComponentImageInit = componentImage;
+    summary.ComponentImageFinished = componentImage;
     %imviewer(componentImage) % Todo: return this as part of summary
 
+    if isempty(S) || isempty(allPixelIndices)
+        % Nothing to pick from; the merge and tagging below need rois.
+        warning('on', 'stats:linkage:NonMonotonicTree')
+        if nargout == 1
+            clear summary
+        end
+        return
+    end
+
     mask = false(imageSize);
 
     nRois = 0;
@@ -83,7 +94,12 @@
 
         % Find peak in the summed component image. This will be the most
         % active roi among candidates
-        [~, peakInd] = max(componentImage(:));
+        [peakValue, peakInd] = max(componentImage(:));
+        if peakValue <= 0
+            % Every remaining pixel has been cleared; without this guard
+            % the loop would spin on an empty peak forever.
+            break
+        end
 
         % Find all components that contain this peak
         containsPeak = regionInd(allPixelIndices==peakInd);

code/+nansen/+module/+twophoton/+integration/+caiman/deconvolve.m

@@ -30,11 +30,7 @@
     % todo...
     prevstr = [];
     starttime = tic;
-    if isempty(gcp('nocreate'))
-        dispProgress = true;
-    else
-        dispProgress = false;
-    end
+    dispProgress = ~isParallelPoolRunning();
 
     if nRois < 5
         dispProgress = false;
@@ -141,3 +137,11 @@
 fprintf(msg)
 end
 end
+
+function tf = isParallelPoolRunning()
+%isParallelPoolRunning True if a parallel pool is open
+%
+%   gcp belongs to the Parallel Computing Toolbox, which the module lists
+%   as optional, so it must not be called when the toolbox is absent.
+    tf = exist('gcp', 'file') == 2 && ~isempty(gcp('nocreate'));
+end

code/+nansen/+module/+twophoton/+integration/+extract/@Processor/Processor.m

@@ -109,20 +109,6 @@ function onCompletion(obj)
                     obj.mergeResults()
                 end
 
-                % % % Todo: The following should be a separate session method...
-                % % if isempty(obj.OriginalStack)
-                % %     sourceStack = obj.SourceStack;
-                % % else
-                % %     sourceStack = obj.OriginalStack;
-                % % end
-                % %
-                % % % Get temporal segments
-                % % tExtracor = nansen.module.twophoton.integration.extract.ProcessorT(...
-                % %     sourceStack, obj.Options, obj.MergedResults);
-                % % tExtracor.Options.Run.numFramesPerPart = 2000;
-                % % tExtracor.DataIoModel = obj.DataIoModel;
-                % %
-                % % tExtracor.runMethod()
             end
 
             %obj.createRoiClassificationData()

code/+nansen/+module/+twophoton/+integration/+suite2p/Suite2pCorrected.m

@@ -13,8 +13,7 @@
     methods % Structors
 
         function obj = Suite2pCorrected(filePath, varargin)
-            import('nansen.module.twophoton.integration.suite2p.Suite2pCorrected')
-            filePath = Suite2pCorrected.lookForMultipartFiles(filePath);
+            filePath = nansen.module.twophoton.integration.suite2p.Suite2pCorrected.lookForMultipartFiles(filePath);
 
             obj@nansen.stack.virtual.TiffMultiPart(filePath, varargin{:})
         end

code/+nansen/+module/+twophoton/+io/+mscan/MDF.m

@@ -122,36 +122,6 @@ function writeFrames(obj, frameIndex, data)
 
     methods (Access = private)
 
-        function resolveDataSizeAndDimensionArrangement(obj, stackSize)
-            % Todo: Remove when this is merged into main branch from
-            % dev-imagestack-5d
-
-            % Find singleton dimensions.
-            isSingleton = stackSize == 1;
-
-            % Get arrangement of dimensions of data
-            try
-                % Note subclasses might implement this as a constant
-                % property. If it's not implemented, use default DDA
-                dataDimensionArrangement = obj.DATA_DIMENSION_ARRANGEMENT;
-            catch
-                dataDimensionArrangement = obj.DEFAULT_DIMENSION_ARRANGEMENT;
-            end
-
-            % Get order of dimensions of data
-            [~, ~, dimensionOrder] = intersect( obj.DEFAULT_DIMENSION_ARRANGEMENT, ...
-                dataDimensionArrangement, 'stable' );
-
-            % Rearrange beased on dimension order
-            isSingleton_(dimensionOrder) = isSingleton;
-            dataSize(dimensionOrder) = stackSize;
-
-            % Assign size and dimension arrangement for data excluding
-            % singleton dimension.
-            obj.DataSize = dataSize(~isSingleton_);
-            obj.DataDimensionArrangement = dataDimensionArrangement(~isSingleton_);
-        end
-
         function mdfParams = getScanParameters(obj)
             import nansen.module.twophoton.io.mscan.getScanParameters
             mdfParams = getScanParameters(obj.mfile);

code/+nansen/+module/+twophoton/+io/+sciscan/SciScanRaw.m

@@ -431,7 +431,7 @@ function getChannelColors(obj)
         inistring = fileread(iniFilePath);
 
         sciscanVars = {'external.start.trigger.enable', 'aocard.model'};
-        isValid = contains(inistring, sciscanVars);
+        isValid = all(contains(inistring, sciscanVars));
     end
 end
 

code/+nansen/+module/+twophoton/+io/+thorlabs/ThorLabsTiffs.m

@@ -43,9 +43,7 @@
 methods % Structors
 
     function obj = ThorLabsTiffs(filePath, varargin)
-        import('nansen.module.twophoton.io.thorlabs.ThorLabsTiffs')
-
-        filePath = ThorLabsTiffs.lookForMultipartFiles(filePath);
+        filePath = nansen.module.twophoton.io.thorlabs.ThorLabsTiffs.lookForMultipartFiles(filePath);
 
         obj@nansen.stack.data.VirtualArray(filePath, varargin{:})
     end

code/+nansen/+module/+twophoton/+roisignals/+process/+dff/dffChenEtAl2013.m

@@ -16,8 +16,8 @@
         fPil = squeeze( mean(signalArray(:, 2:end, :), 2) );
     end
 
-    npil_true0 = prctile(fPil, 20);
+    npil_true0 = prctile(fPil, 20, 1);
     f_true = fRoi - (0.7*fPil) + npil_true0;
-    f_true0 = prctile(f_true, 20);
+    f_true0 = prctile(f_true, 20, 1);
     dff_true = (f_true - f_true0) ./ f_true0;
 end

code/+nansen/+module/+twophoton/+roisignals/+process/+dff/dffRoiMinusDffNpil.m

@@ -51,6 +51,13 @@
         if ~isa(deltaDff, 'double') || ~isa(dffNpil, 'double')
             deltaDff = double(deltaDff); dffNpil = double(dffNpil);
         end
+        if exist('sgolayfilt', 'file') ~= 2
+            error('NANSEN:TwoPhoton:SignalProcessingToolboxRequired', ...
+                ['The dF/F method "dffRoiMinusDffNpil" uses sgolayfilt from the ', ...
+                 'Signal Processing Toolbox, which is not installed. Install the ', ...
+                 'toolbox, set smoothSignals to true, or choose another dffFcn ', ...
+                 'such as dffClassic.'])
+        end
         baselineCorrection = sgolayfilt(deltaDff, 3, 11);
         dff = dffRoi - sgolayfilt(dffNpil, 3, 11) + baselineCorrection;
     end

code/+nansen/+module/+twophoton/+roisignals/deconvolveDff.m

@@ -1,4 +1,20 @@
 function [dec, den, opt] = deconvolveDff(dff, varargin)
+%deconvolveDff Deconvolve dF/F signals with the selected method
+%
+%   [dec, den, opt] = deconvolveDff(dff) deconvolves dff, a numSamples x
+%   numRois array as returned by computeDff, and returns the deconvolved
+%   (dec) and denoised (den) signals in the same shape, plus the options
+%   used (opt). A vector is treated as a single ROI in either orientation.
+%
+%   [dec, den, opt] = deconvolveDff(dff, Name, Value, ...) or
+%   deconvolveDff(dff, options) sets the deconvolution parameters, see
+%   getDeconvolutionParameters. The deconvolutionMethod parameter selects
+%   a function from the process.deconvolve package by name.
+%
+%   P = deconvolveDff() returns the default parameters.
+%
+%   See also nansen.module.twophoton.roisignals.computeDff,
+%   nansen.module.twophoton.roisignals.getDeconvolutionParameters
 
     [P, V] = nansen.module.twophoton.roisignals.getDeconvolutionParameters();
     P.deconvolutionMethod = 'caiman';
@@ -12,17 +28,19 @@
     deconvPackage = 'nansen.module.twophoton.roisignals.process.deconvolve';
     deconvFunction = str2func( strjoin({deconvPackage, params.deconvolutionMethod}, '.') );
 
-    % dff must be nRois x nSamples for deconvolution method...
-    wasTransposed = false;
-    if size(dff, 1) > size(dff, 2)
-        dff = transpose(dff);
-        wasTransposed = true;
+    % The deconvolution methods want numRois x numSamples. The input
+    % contract is numSamples x numRois (computeDff's output); a vector is
+    % one ROI regardless of orientation. Shape is not guessed from which
+    % dimension is larger, so recordings with more ROIs than samples are
+    % deconvolved correctly.
+    if isvector(dff)
+        dffForMethod = reshape(dff, 1, []);
+    else
+        dffForMethod = transpose(dff);
     end
 
-    [dec, den, opt] = deconvFunction(dff, params);
+    [dec, den, opt] = deconvFunction(dffForMethod, params);
 
-    if wasTransposed
-        dec = transpose(dec);
-        den = transpose(den);
-    end
+    dec = reshape(transpose(dec), size(dff));
+    den = reshape(transpose(den), size(dff));
 end

code/+nansen/+module/+twophoton/+roisignals/extractF.m

@@ -3,7 +3,8 @@
 %
 %   signals = extractF(imageStack, roiArray) extracts signals for rois in
 %   roiArray from image frames in imageData using default settings.
-%   roiArray must be a nansen.RoiArray object and imageData can be a
+%   roiArray must be an array of RoI objects (or a roimanager.roiGroup)
+%   and imageData can be a
 %   numeric array of three or more dimensions or a nansen.ImageStack
 %   object.
 %