diff --git a/.github/badges/code_issues.svg b/.github/badges/code_issues.svg
index edc4000eb..e0d2558a1 100644
--- a/.github/badges/code_issues.svg
+++ b/.github/badges/code_issues.svg
@@ -1 +1 @@
-
+
diff --git a/.github/badges/tests.svg b/.github/badges/tests.svg
index fd00e9b16..783fa2a5e 100644
--- a/.github/badges/tests.svg
+++ b/.github/badges/tests.svg
@@ -1 +1 @@
-
+
diff --git a/code/+nansen/+processing/+roi/createMasks.m b/code/+nansen/+processing/+roi/createMasks.m
index 8b50e4fa0..7be319e74 100644
--- a/code/+nansen/+processing/+roi/createMasks.m
+++ b/code/+nansen/+processing/+roi/createMasks.m
@@ -50,7 +50,12 @@
return
end
- neuropilMaskArray = estimateNeuropilMasks(roiMaskArrayOrig, params);
+ neuropilMaskArray = estimateNeuropilMasks(roiMaskArrayOrig, ...
+ roiInd=params.roiInd, ...
+ imageMask=params.imageMask, ...
+ excludeRoiFromNeuropil=params.excludeRoiFromNeuropil, ...
+ neuropilExpansionFactor=params.neuropilExpansionFactor, ...
+ cellNeuropilSeparation=params.cellNeuropilSeparation);
if params.numNeuropilSlices > 1
neuropilMasks = splitNeuropilMasks(neuropilMaskArray, ...
diff --git a/code/+nansen/+processing/+roi/estimateNeuropilMasks.m b/code/+nansen/+processing/+roi/estimateNeuropilMasks.m
index 3dc2d0bb8..a8858191f 100644
--- a/code/+nansen/+processing/+roi/estimateNeuropilMasks.m
+++ b/code/+nansen/+processing/+roi/estimateNeuropilMasks.m
@@ -1,36 +1,105 @@
-function neuropilMasks = estimateNeuropilMasks(roiMasksIn, varargin)
-%estimateNeuropilMasks Estimate neuropil masks for RoIs.
+function neuropilMasks = estimateNeuropilMasks(roiMasks, options)
+%estimateNeuropilMasks - Create a mask of the neuropil around each roi
+% neuropilMasks = estimateNeuropilMasks(roiMasks) creates one neuropil
+% mask for each roi mask. roiMasks is a logical array of size
+% imageHeight-by-imageWidth-by-numRois and neuropilMasks is a logical
+% array of the same size. A neuropil mask is grown from its roi mask in
+% steps of one pixel. The first step adds the four edge neighbors of
+% each mask pixel, the second step adds the four corner neighbors, and
+% the following steps alternate between the two. Growing stops when
+% the mask holds at least neuropilExpansionFactor times as many pixels
+% as the roi, or after 101 steps. The pixels of a roi are never part of
+% its own neuropil mask. Pixels that imageMask or excludeRoiFromNeuropil
+% leave out of the neuropil masks are not counted.
%
+% neuropilMasks = estimateNeuropilMasks(roiMasks,Name=VALUE) specifies
+% options using one or more name-value arguments. An option that is not
+% given takes its default value from
+% nansen.twophoton.roisignals.extract.getDefaultParameters:
+% roiInd - Indices of the rois to create neuropil
+% masks for, or "all" (default). The
+% third dimension of neuropilMasks
+% follows the order of roiInd.
+% imageMask - Logical matrix of size
+% imageHeight-by-imageWidth. Pixels that
+% are false are left out of all neuropil
+% masks. The default is empty, which
+% leaves out no pixels.
+% excludeRoiFromNeuropil - True (default) leaves the pixels of
+% all rois in roiMasks out of all
+% neuropil masks.
+% neuropilExpansionFactor - Number of pixels in a neuropil mask
+% relative to the number of pixels in its
+% roi. The default is 4.
+% cellNeuropilSeparation - Number of one-pixel steps that the
+% rois are grown by before their pixels
+% are left out of the neuropil masks. It
+% has no effect if excludeRoiFromNeuropil
+% is false. The default is 0.
+%
+% Example: Neuropil mask around one of two neighboring rois
+% roiMasks = false(64,64,2);
+% roiMasks(20:30,20:30,1) = true;
+% roiMasks(20:30,34:44,2) = true;
+% neuropilMasks = nansen.processing.roi.estimateNeuropilMasks( ...
+% roiMasks,roiInd=1);
+% imagesc(neuropilMasks + 2*any(roiMasks,3))
+%
+% See also nansen.processing.roi.createMasks,
+% nansen.processing.roi.prepareRoiMasks,
+% nansen.twophoton.roisignals.extractF,
+% nansen.twophoton.roisignals.extract.getDefaultParameters
+
% TODO
- % [ ] Instead of using max count of 100 to stop, use a parameter
- % called max neuropil radius....
- % [ ] Add doc
- % [ ] test the max radius
- % [ ] add warning if mask can't be created (i.e not enough neuropil available)
+ % [ ] Stop growing at the maxNeuropilRadius parameter instead of
+ % after a fixed number of steps.
+ % [ ] Add warning if mask can't be created (i.e not enough neuropil
+ % available)
+
+ arguments
+ roiMasks logical
+ options.roiInd {mustBeRoiIndices}
+ options.imageMask (:,:) logical
+ options.excludeRoiFromNeuropil (1,1) logical
+ options.neuropilExpansionFactor (1,1) {mustBeNumeric, mustBeReal, mustBePositive}
+ options.cellNeuropilSeparation (1,1) {mustBeNumeric, mustBeReal, mustBeInteger, mustBeNonnegative}
+ end
+
+ [height, width, numRois] = size(roiMasks);
- % Get default parameters and assertion functions.
- [P, V] = nansen.twophoton.roisignals.extract.getDefaultParameters();
+ params = addDefaultValues(options);
- % Parse potential parameters from input arguments
- params = utility.parsenvpairs(P, V, varargin{:});
+ if isIndexOfAllRois(params.roiInd)
+ params.roiInd = 1:numRois;
+ end
+
+ hasImageMask = ~isempty(params.imageMask);
+ if hasImageMask && ~isequal(size(params.imageMask), [height, width])
+ error("NANSEN:RoiSignals:ImageMaskSizeMismatch", ...
+ "The size of imageMask is %s, but the size of the roi masks is %s. " + ...
+ "Give an imageMask with one element for each pixel of a roi mask.", ...
+ mat2str(size(params.imageMask)), mat2str([height, width]))
+ end
% Create mask of pixel to ignore when creating neuropil masks
if params.excludeRoiFromNeuropil
% Create an array with all rois overlapping.
- footprintAllRois = sum(roiMasksIn, 3) >= 1;
+ footprintAllRois = sum(roiMasks, 3) >= 1;
- % Use binary erosion to remove pixels in the vicinity of main rois
+ % Use binary dilation to add pixels in the vicinity of main rois
for i = 1:params.cellNeuropilSeparation
footprintAllRois = imdilateoddeven(footprintAllRois, i);
end
- if ~isempty(params.imageMask)
- pixelsToIgnore = footprintAllRois & ~params.imageMask;
+ % Ignore a pixel if it belongs to a roi or if it is excluded by
+ % the image mask
+ if hasImageMask
+ pixelsToIgnore = footprintAllRois | ~params.imageMask;
else
pixelsToIgnore = footprintAllRois;
end
else
- if ~isempty(params.imageMask)
+ if hasImageMask
pixelsToIgnore = ~params.imageMask;
else
pixelsToIgnore = [];
@@ -38,14 +107,13 @@
end
% Allocate array for neuropil masks.
- [height, width, ~] = size(roiMasksIn);
neuropilMasks = false( height, width, numel(params.roiInd) );
% Loop through rois
for n = 1:numel(params.roiInd)
currentRoiInd = params.roiInd(n);
- currentRoiMask = roiMasksIn(:, :, currentRoiInd);
+ currentRoiMask = roiMasks(:, :, currentRoiInd);
thisMask = currentRoiMask;
origArea = sum(currentRoiMask(:)); % original area
@@ -93,6 +161,37 @@
end
end
+function params = addDefaultValues(options)
+%addDefaultValues - Add the default value of each option that is not given
+ defaults = nansen.twophoton.roisignals.extract.getDefaultParameters();
+
+ optionNames = ["roiInd", "imageMask", "excludeRoiFromNeuropil", ...
+ "neuropilExpansionFactor", "cellNeuropilSeparation"];
+
+ params = options;
+ for name = optionNames
+ if ~isfield(params, name)
+ params.(name) = defaults.(name);
+ end
+ end
+end
+
+function tf = isIndexOfAllRois(roiInd)
+%isIndexOfAllRois - Check if roiInd is the text "all"
+ tf = (ischar(roiInd) || isstring(roiInd)) && isequal(string(roiInd), "all");
+end
+
+function mustBeRoiIndices(value)
+%mustBeRoiIndices - Validate that value is "all" or a vector of indices
+ isIndexVector = isnumeric(value) && isvector(value) && isreal(value) ...
+ && all(value == round(value)) && all(value >= 1);
+
+ if ~(isIndexOfAllRois(value) || isIndexVector)
+ error("NANSEN:RoiSignals:InvalidRoiIndices", ...
+ "Value must be ""all"" or a vector of positive integers.")
+ end
+end
+
function mask = imdilateoddeven(mask, count)
%imdilateoddeven imdilate using different nhood on odd and even iterations
%
diff --git a/code/+nansen/+twophoton/+roisignals/extractF.m b/code/+nansen/+twophoton/+roisignals/extractF.m
index 86a327359..b37a7289a 100644
--- a/code/+nansen/+twophoton/+roisignals/extractF.m
+++ b/code/+nansen/+twophoton/+roisignals/extractF.m
@@ -1,33 +1,105 @@
-function [signalArray, P] = extractF(imageData, roiData, varargin)
-%extractF Extract RoI signals from an image stack in chunkwise manner
+function [signalArray, params] = extractF(imageData, roiData, optionsStruct, options)
+%extractF - Extract roi signals from an image stack
+% signalArray = extractF(imageData,roiData) computes the mean pixel
+% value of each roi in roiData for every frame of imageData. imageData
+% is a nansen.stack.ImageStack or a numeric array of size
+% imageHeight-by-imageWidth-by-numFrames. roiData is an array of RoI
+% objects or a roimanager.roiGroup, and the imagesize property of the
+% rois must equal [imageHeight, imageWidth]. signalArray is of size
+% numFrames-by-numSubregions-by-numRois. signalArray(:,1,:) holds the
+% signals of the rois and signalArray(:,2:end,:) holds the signals of
+% the neuropil masks that are created around each roi. Frames are read
+% and processed in chunks if imageData is a virtual stack or if there
+% are 100 rois or more. The chooseChunkLength method of
+% nansen.stack.ImageStack sets the number of frames in a chunk.
%
-% 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
-% numeric array of three or more dimensions or a nansen.ImageStack
-% object.
+% signalArray = extractF(imageData,roiData,optionsStruct) specifies
+% options as fields of the scalar struct optionsStruct. The field names
+% are the names of the name-value arguments below. Fields with other
+% names are ignored, so that optionsStruct can be the struct returned
+% by nansen.twophoton.roisignals.extract.getDefaultParameters.
%
-% signals = extractF(imageStack, roiArray, params) extracts signals
-% using options specified in params. Params gan be given as a
-% struct or as a cell array of name-value pairs.
+% signalArray = extractF(...,Name=VALUE) specifies options using one or
+% more name-value arguments. An option given by name replaces the field
+% of the same name in optionsStruct. An option that is not given takes
+% its default value from
+% nansen.twophoton.roisignals.extract.getDefaultParameters:
+% roiInd - Indices of the rois to extract signals
+% for, or "all" (default).
+% imageMask - Logical matrix of size
+% imageHeight-by-imageWidth. Pixels that
+% are false are left out of the neuropil
+% masks. The roi masks are not affected.
+% The default includes all pixels.
+% pixelComputationMethod - Statistic that is computed over the
+% pixels of a mask, "mean" (default) or
+% "median". "median" is always computed
+% by serialExtract.
+% excludeRoiOverlaps - True (default) leaves pixels that
+% belong to more than one roi out of the
+% roi masks.
+% createNeuropilMask - True (default) creates a neuropil mask
+% around each roi and extracts its
+% signal.
+% excludeRoiFromNeuropil - True (default) leaves the pixels of
+% all rois out of the neuropil masks.
+% neuropilExpansionFactor - Number of pixels in a neuropil mask
+% relative to the number of pixels in
+% its roi. The default is 4.
+% cellNeuropilSeparation - Number of one-pixel steps that the
+% rois are grown by before their pixels
+% are left out of the neuropil masks.
+% The default is 0.
+% numNeuropilSlices - Number of slices that each neuropil
+% mask is split into. Each slice is one
+% subregion in signalArray. The default
+% is 1.
+% maxNeuropilRadius - Not used.
+% extractFcn - Handle of the function that extracts
+% the signals from a chunk of frames.
+% The default is serialExtract for fewer
+% than 200 rois and batchExtract for 200
+% rois or more. Both are in the package
+% nansen.twophoton.roisignals.extract.
+% roiMaskFormat - Format of the prepared roi masks,
+% "struct" for serialExtract and
+% "sparse" for batchExtract. It is set
+% to match extractFcn.
+% signalDataType - Class of signalArray, "single"
+% (default) or "double".
+% channelNum - Channels to extract signals from if
+% imageData has more than one current
+% channel. The signals of the channels
+% are concatenated along the third
+% dimension of signalArray. If roiData
+% is an array of roimanager.roiGroup
+% objects, roiData(i) is used for
+% channel i. The default is 1.
+% verbose - True prints a message for each chunk
+% of frames and the elapsed time. The
+% default is false.
+% showTimer - True prints the elapsed time. The
+% default is false.
%
-% Parameters:
+% [signalArray,PARAMS] = extractF(...) also returns the options that
+% were used, as a scalar struct with one field for each name-value
+% argument. PARAMS.roiInd holds the indices of the rois,
+% PARAMS.imageMask the image mask that was applied, and
+% PARAMS.extractFcn and PARAMS.roiMaskFormat the extraction function
+% and the mask format that were selected. If imageData has more than
+% one current channel, these four fields hold the values that were
+% given, because they are resolved separately for each channel.
%
-% roiInd : integer vector in range [1, numRois]
-% List of roi indices (use for extraction of signals from a subset
-% of rois). Default : extract signals from all rois.
-% imageMask : logical matrix (imageHeight x imageWidth)
-% Use for excluding regions of image. Include pixels that are true
-% and exclude pixels that are false. Default : include all pixels
-% excludeRoiOverlaps : logical scalar
-% Exclude pixels where rois are overlapping. Default = true
-% createNeuropilMask : logical scalar
-% Create mask (and extract signals) for surrounding neuropil regions
-% Default = true
-% excludeRoiFromNeuropil : logical scalar
-% Exclude rois from neuropil regions. Default = true
+% Example: Extract the signal of one roi from a numeric array
+% imageData = rand(64,64,100);
+% roiArray = RoI("Circle",[32,32,5],[64,64]);
+% signalArray = nansen.twophoton.roisignals.extractF( ...
+% imageData,roiArray,createNeuropilMask=false);
+% plot(signalArray)
%
-% See also nansen.twophoton.roisignals.extract.getDefaultParameters
+% See also nansen.twophoton.roisignals.computeDff,
+% nansen.twophoton.roisignals.extract.getDefaultParameters,
+% nansen.processing.roi.prepareRoiMasks, nansen.stack.ImageStack, RoI
% Todo:
% [x] Implement a standard set of options across signal extraction
@@ -37,9 +109,9 @@
% versions of signals????
%
% [x] Make sure it works to extract from subset of rois....
- % [ ] Test that it works as expected. I.e do we get signals from
- % requested rois, and are masks created according to options. Use
- % imviewer/roimanager gui for this...
+ % [x] Test that it works as expected. I.e do we get signals from
+ % requested rois, and are masks created according to options. See
+ % nansen.unittest.roisignals.ExtractFTest
%
% [ ] Extraction function should depend on number of rois. If nRois <
% 140 should use looping over rois and array cropping, otherwise
@@ -50,16 +122,44 @@
%
% [ ] Support for multiple planes.
- % Get default parameters and assertion functions.
+ arguments
+ imageData
+ roiData
+ optionsStruct struct {mustBeScalarOrEmpty} = struct.empty
+ options.roiInd {mustBeRoiIndices}
+ options.imageMask (:,:) logical
+ options.pixelComputationMethod {mustBeTextScalar, ...
+ mustBeMember(options.pixelComputationMethod, ["mean", "median"])}
+ options.excludeRoiOverlaps (1,1) logical
+ options.createNeuropilMask (1,1) logical
+ options.excludeRoiFromNeuropil (1,1) logical
+ options.neuropilExpansionFactor (1,1) {mustBeNumeric, mustBeReal, mustBePositive}
+ options.cellNeuropilSeparation (1,1) {mustBeNumeric, mustBeReal, mustBeInteger, mustBeNonnegative}
+ options.numNeuropilSlices (1,1) {mustBeNumeric, mustBeReal, mustBeInteger, mustBeNonnegative}
+ options.maxNeuropilRadius (1,1) {mustBeNumeric, mustBeReal, mustBeInteger, mustBeNonnegative}
+ options.roiMaskFormat {mustBeTextScalar, ...
+ mustBeMember(options.roiMaskFormat, ["sparse", "struct"])}
+ options.extractFcn {mustBeFunctionHandleOrEmpty}
+ options.showTimer (1,1) logical
+ options.verbose (1,1) logical
+ options.signalDataType {mustBeTextScalar, ...
+ mustBeMember(options.signalDataType, ["single", "double"])}
+ options.channelNum {mustBeNumeric, mustBeReal, mustBeInteger, mustBePositive, mustBeNonempty}
+ end
- [P, V] = nansen.twophoton.roisignals.extract.getDefaultParameters();
- P.showTimer = false; V.showTimer = @(x) assert(islogical(x), 'Value must be logical');
- P.verbose = false; V.verbose = @(x) assert(islogical(x), 'Value must be logical');
- P.signalDataType = 'single'; V.signalDataType = @(x) assert(any(strcmp(x, {'single', 'double'})), 'Value must be ''single'' or ''double''');
- P.channelNum = 1; V.channelNum = @(x) not(isempty(x));
+ % Options in a struct are passed on as name-value arguments, so that
+ % they are validated in the same way as options given by name. No
+ % option has a default value in the arguments block, because a default
+ % value could not be told apart from a value given by name, and would
+ % replace the value in the struct.
+ if ~isempty(optionsStruct)
+ nvPairs = combineOptions(optionsStruct, options); % Local function
+ [signalArray, params] = nansen.twophoton.roisignals.extractF( ...
+ imageData, roiData, nvPairs{:});
+ return
+ end
- % Parse potential parameters from input arguments
- params = utility.parsenvpairs(P, V, varargin{:});
+ params = addDefaultValues(options); % Local function
% Validate roidata
if isa(roiData, 'roimanager.roiGroup')
@@ -79,13 +179,14 @@
% If multiple channels are selected, extract from each individually
numChannels = numel( imageStack.CurrentChannel );
if numChannels > 1
+ nvPairs = namedargs2cell(options);
if ~isempty(params.channelNum)
currentChannel = imageStack.CurrentChannel;
imageStack.CurrentChannel = params.channelNum;
- signalArray = extractFromMultipleChannels(imageStack, roiData, varargin{:});
+ signalArray = extractFromMultipleChannels(imageStack, roiData, nvPairs{:});
imageStack.CurrentChannel = currentChannel;
else
- signalArray = extractFromMultipleChannels(imageStack, roiData, varargin{:});
+ signalArray = extractFromMultipleChannels(imageStack, roiData, nvPairs{:});
end
return
end
@@ -121,7 +222,6 @@
[elapsedTimeLoad, elapsedTimeExtract] = deal( 0 );
signalExtractionFcn = params.extractFcn;
- method = params.pixelComputationMethod;
if params.verbose
fprintf('Image stack is split in %d parts for signal extraction\n', numParts)
@@ -151,6 +251,69 @@
end
end
+function params = addDefaultValues(options)
+%addDefaultValues - Add the default value of each option that is not given
+ params = nansen.twophoton.roisignals.extract.getDefaultParameters();
+
+ % Remove the fields that configure the options editor. The name of
+ % such a field is the name of a parameter followed by an underscore.
+ names = string(fieldnames(params));
+ isConfigField = endsWith(names, "_") & ...
+ ismember(extractBefore(names, strlength(names)), names);
+ params = rmfield(params, names(isConfigField));
+
+ % Options of extractF that are not parameters of the extraction
+ % functions
+ params.showTimer = false;
+ params.verbose = false;
+ params.signalDataType = "single";
+ params.channelNum = 1;
+
+ names = fieldnames(options);
+ for i = 1:numel(names)
+ params.(names{i}) = options.(names{i});
+ end
+end
+
+function nvPairs = combineOptions(optionsStruct, options)
+%combineOptions - Combine an options struct with name-value arguments
+% Fields of optionsStruct that are not options of extractF are left out.
+% A field of options replaces the field of the same name in
+% optionsStruct.
+ optionNames = fieldnames(addDefaultValues(struct()));
+
+ fieldNames = fieldnames(optionsStruct);
+ isOption = ismember(fieldNames, optionNames);
+ combined = rmfield(optionsStruct, fieldNames(~isOption));
+
+ names = fieldnames(options);
+ for i = 1:numel(names)
+ combined.(names{i}) = options.(names{i});
+ end
+
+ nvPairs = namedargs2cell(combined);
+end
+
+function mustBeRoiIndices(value)
+%mustBeRoiIndices - Validate that value is "all" or a vector of indices
+ isAll = (ischar(value) || isstring(value)) && isequal(string(value), "all");
+ isIndexVector = isnumeric(value) && isvector(value) && isreal(value) ...
+ && all(value == round(value)) && all(value >= 1);
+
+ if ~(isAll || isIndexVector)
+ error("NANSEN:RoiSignals:InvalidRoiIndices", ...
+ "Value must be ""all"" or a vector of positive integers.")
+ end
+end
+
+function mustBeFunctionHandleOrEmpty(value)
+%mustBeFunctionHandleOrEmpty - Validate that value is a function handle or empty
+ if ~(isempty(value) || isa(value, "function_handle"))
+ error("NANSEN:RoiSignals:InvalidExtractFcn", ...
+ "Value must be a function handle or empty.")
+ end
+end
+
function validateInputDimensions(imageStack, roiArray)
%validateInputDimensions Check that ImageStack and roiArray has matching
% dimensions.
diff --git a/code/modules/+nansen/+module/+ophys/+twophoton/+sessionmethod/+process/+signalExtraction/extractSignals.m b/code/modules/+nansen/+module/+ophys/+twophoton/+sessionmethod/+process/+signalExtraction/extractSignals.m
index 6b9b7e6df..df6002255 100644
--- a/code/modules/+nansen/+module/+ophys/+twophoton/+sessionmethod/+process/+signalExtraction/extractSignals.m
+++ b/code/modules/+nansen/+module/+ophys/+twophoton/+sessionmethod/+process/+signalExtraction/extractSignals.m
@@ -67,7 +67,7 @@ function runMethod(obj)
roiArray = sessionData.RoiArray;
extractF = @nansen.twophoton.roisignals.extractF;
- [signalArray, P] = extractF(imageStack, roiArray, 'verbose', true, obj.Options);
+ [signalArray, P] = extractF(imageStack, roiArray, obj.Options, verbose=true);
% Todo: Save results...
obj.saveData('RoiSignals_MeanF', squeeze(signalArray(:, 1, :)) )
diff --git a/tests/+nansen/+unittest/+roisignals/EstimateNeuropilMasksTest.m b/tests/+nansen/+unittest/+roisignals/EstimateNeuropilMasksTest.m
new file mode 100644
index 000000000..da67f516b
--- /dev/null
+++ b/tests/+nansen/+unittest/+roisignals/EstimateNeuropilMasksTest.m
@@ -0,0 +1,157 @@
+classdef EstimateNeuropilMasksTest < matlab.unittest.TestCase
+ %EstimateNeuropilMasksTest Unit tests for estimateNeuropilMasks.
+ %
+ % Tests nansen.processing.roi.estimateNeuropilMasks on three square
+ % roi masks in a 64-by-64 image. Roi 1 and roi 2 are 3 pixels apart,
+ % so that roi 2 lies where the neuropil mask of roi 1 grows. Roi 3
+ % has no neighbors.
+ %
+ % testMasksLeaveOutRoiPixels with the allTrue parameter is a
+ % regression test: an image mask that includes all pixels used to
+ % switch off excludeRoiFromNeuropil, and extractF always passes such
+ % a mask.
+ %
+ % Run tests:
+ % runtests('nansen.unittest.roisignals.EstimateNeuropilMasksTest')
+
+ properties (Constant, Access = private)
+ ImageSize = [64, 64]
+ end
+
+ properties (TestParameter)
+ % Image masks that leave out no pixels
+ fullImageMask = struct( ...
+ 'empty', {logical.empty}, ...
+ 'allTrue', {true(64, 64)})
+ end
+
+ methods (Test)
+
+ function testMasksLeaveOutRoiPixels(testCase, fullImageMask)
+ roiMasks = createRoiMasks(testCase.ImageSize);
+
+ neuropilMasks = nansen.processing.roi.estimateNeuropilMasks( ...
+ roiMasks, imageMask=fullImageMask);
+
+ isRoiPixel = any(roiMasks, 3);
+ testCase.verifyFalse(any(neuropilMasks & isRoiPixel, "all"))
+ end
+
+ function testNeighborRoiIsKeptWhenRoisAreNotExcluded(testCase)
+ roiMasks = createRoiMasks(testCase.ImageSize);
+
+ neuropilMasks = nansen.processing.roi.estimateNeuropilMasks( ...
+ roiMasks, excludeRoiFromNeuropil=false);
+
+ maskOfRoi1 = neuropilMasks(:, :, 1);
+ testCase.verifyTrue(any(maskOfRoi1 & roiMasks(:, :, 2), "all"), ...
+ "The neuropil mask of roi 1 should reach into roi 2.")
+ testCase.verifyFalse(any(maskOfRoi1 & roiMasks(:, :, 1), "all"), ...
+ "A neuropil mask should never hold pixels of its own roi.")
+ end
+
+ function testMasksLeaveOutPixelsOutsideImageMask(testCase)
+ roiMasks = createRoiMasks(testCase.ImageSize);
+ imageMask = true(testCase.ImageSize);
+ imageMask(:, 1:19) = false; % The columns to the left of roi 1
+
+ neuropilMasks = nansen.processing.roi.estimateNeuropilMasks( ...
+ roiMasks, imageMask=imageMask);
+
+ testCase.verifyFalse(any(neuropilMasks & ~imageMask, "all"))
+ testCase.verifyTrue(any(neuropilMasks(:, :, 1), "all"))
+ end
+
+ function testImageMaskAppliesWhenRoisAreNotExcluded(testCase)
+ roiMasks = createRoiMasks(testCase.ImageSize);
+ imageMask = true(testCase.ImageSize);
+ imageMask(:, 1:19) = false;
+
+ neuropilMasks = nansen.processing.roi.estimateNeuropilMasks( ...
+ roiMasks, imageMask=imageMask, excludeRoiFromNeuropil=false);
+
+ testCase.verifyFalse(any(neuropilMasks & ~imageMask, "all"))
+ end
+
+ function testMaskAreaReachesExpansionFactor(testCase)
+ roiMasks = createRoiMasks(testCase.ImageSize);
+ expansionFactor = 2;
+
+ neuropilMasks = nansen.processing.roi.estimateNeuropilMasks( ...
+ roiMasks, neuropilExpansionFactor=expansionFactor);
+
+ roiArea = squeeze(sum(roiMasks, [1, 2]));
+ neuropilArea = squeeze(sum(neuropilMasks, [1, 2]));
+ testCase.verifyGreaterThanOrEqual(neuropilArea, expansionFactor * roiArea)
+ end
+
+ function testLargerExpansionFactorGivesLargerMasks(testCase)
+ roiMasks = createRoiMasks(testCase.ImageSize);
+
+ smallMasks = nansen.processing.roi.estimateNeuropilMasks( ...
+ roiMasks, neuropilExpansionFactor=2);
+ largeMasks = nansen.processing.roi.estimateNeuropilMasks( ...
+ roiMasks, neuropilExpansionFactor=4);
+
+ smallArea = squeeze(sum(smallMasks, [1, 2]));
+ largeArea = squeeze(sum(largeMasks, [1, 2]));
+ testCase.verifyGreaterThan(largeArea, smallArea)
+ end
+
+ function testSeparationKeepsMasksAwayFromRois(testCase)
+ roiMasks = createRoiMasks(testCase.ImageSize);
+
+ neuropilMasks = nansen.processing.roi.estimateNeuropilMasks( ...
+ roiMasks, cellNeuropilSeparation=2);
+
+ % Two steps cover the eight neighbors of every roi pixel
+ isNextToRoi = imdilate(any(roiMasks, 3), ones(3));
+ testCase.verifyFalse(any(neuropilMasks & isNextToRoi, "all"))
+ end
+
+ function testRoiIndSelectsAndOrdersMasks(testCase)
+ roiMasks = createRoiMasks(testCase.ImageSize);
+
+ allMasks = nansen.processing.roi.estimateNeuropilMasks(roiMasks);
+ selectedMasks = nansen.processing.roi.estimateNeuropilMasks( ...
+ roiMasks, roiInd=[3, 1]);
+
+ testCase.verifyEqual(selectedMasks, allMasks(:, :, [3, 1]))
+ end
+
+ function testRoiIndAllEqualsDefault(testCase)
+ roiMasks = createRoiMasks(testCase.ImageSize);
+
+ defaultMasks = nansen.processing.roi.estimateNeuropilMasks(roiMasks);
+ allMasks = nansen.processing.roi.estimateNeuropilMasks( ...
+ roiMasks, roiInd="all");
+
+ testCase.verifySize(allMasks, size(roiMasks))
+ testCase.verifyEqual(allMasks, defaultMasks)
+ end
+
+ function testImageMaskOfWrongSizeThrowsError(testCase)
+ roiMasks = createRoiMasks(testCase.ImageSize);
+
+ testCase.verifyError( ...
+ @() nansen.processing.roi.estimateNeuropilMasks(roiMasks, imageMask=true(10, 10)), ...
+ "NANSEN:RoiSignals:ImageMaskSizeMismatch")
+ end
+
+ function testInvalidRoiIndThrowsError(testCase)
+ roiMasks = createRoiMasks(testCase.ImageSize);
+
+ testCase.verifyError( ...
+ @() nansen.processing.roi.estimateNeuropilMasks(roiMasks, roiInd=0), ...
+ "NANSEN:RoiSignals:InvalidRoiIndices")
+ end
+ end
+end
+
+function roiMasks = createRoiMasks(imageSize)
+%createRoiMasks Two neighboring rois and one roi without neighbors
+ roiMasks = false([imageSize, 3]);
+ roiMasks(20:30, 20:30, 1) = true;
+ roiMasks(20:30, 34:44, 2) = true;
+ roiMasks(45:50, 8:13, 3) = true;
+end
diff --git a/tests/+nansen/+unittest/+roisignals/ExtractFTest.m b/tests/+nansen/+unittest/+roisignals/ExtractFTest.m
new file mode 100644
index 000000000..ef54c3783
--- /dev/null
+++ b/tests/+nansen/+unittest/+roisignals/ExtractFTest.m
@@ -0,0 +1,227 @@
+classdef ExtractFTest < matlab.unittest.TestCase
+ %ExtractFTest Unit tests for nansen.twophoton.roisignals.extractF.
+ %
+ % The extracted signals are compared with the mean pixel values of
+ % hand-made masks. Three square rois lie in a 48-by-64 image: roi 1
+ % and roi 2 overlap in three columns and roi 3 has no neighbors.
+ %
+ % Two tests are regression tests. In
+ % testNeuropilSignalLeavesOutNeighboringRois the rois are brighter
+ % than the background, so the neuropil signal only equals the
+ % background if no neuropil mask holds a pixel of any roi. The
+ % testSecondOutput* tests check that the second output holds the
+ % options that were used and not the default options.
+ %
+ % Run tests:
+ % runtests('nansen.unittest.roisignals.ExtractFTest')
+
+ properties (Constant, Access = private)
+ ImageSize = [48, 64]
+ NumFrames = 25
+ AbsTol = 1e-9 % Signals are requested as double
+ end
+
+ methods (TestMethodSetup)
+ function seedRandomNumberGenerator(testCase)
+ originalState = rng();
+ testCase.addTeardown(@() rng(originalState))
+ rng(1, "twister")
+ end
+ end
+
+ methods (Test)
+
+ function testRoiSignalIsMeanOfNonOverlappingPixels(testCase)
+ [rois, roiMasks] = createRois(testCase.ImageSize);
+ imageData = createRandomStack(testCase.ImageSize, testCase.NumFrames);
+ isInOneRoi = sum(roiMasks, 3) == 1;
+
+ signalArray = nansen.twophoton.roisignals.extractF( ...
+ imageData, rois, signalDataType="double");
+
+ expected = meanInMasks(imageData, roiMasks & isInOneRoi);
+ testCase.verifySize(signalArray, [testCase.NumFrames, 2, 3])
+ testCase.verifyEqual(squeeze(signalArray(:, 1, :)), expected, ...
+ AbsTol=testCase.AbsTol)
+ end
+
+ function testRoiSignalUsesOverlappingPixelsOnRequest(testCase)
+ [rois, roiMasks] = createRois(testCase.ImageSize);
+ imageData = createRandomStack(testCase.ImageSize, testCase.NumFrames);
+
+ signalArray = nansen.twophoton.roisignals.extractF( ...
+ imageData, rois, excludeRoiOverlaps=false, signalDataType="double");
+
+ expected = meanInMasks(imageData, roiMasks);
+ testCase.verifyEqual(squeeze(signalArray(:, 1, :)), expected, ...
+ AbsTol=testCase.AbsTol)
+ end
+
+ function testNeuropilSignalLeavesOutNeighboringRois(testCase)
+ [rois, roiMasks] = createRois(testCase.ImageSize);
+ backgroundValue = 1000;
+ imageData = createTwoLevelStack(any(roiMasks, 3), ...
+ testCase.NumFrames, backgroundValue, 6000);
+
+ signalArray = nansen.twophoton.roisignals.extractF( ...
+ imageData, rois, signalDataType="double");
+
+ expected = backgroundValue * ones(testCase.NumFrames, 3);
+ testCase.verifyEqual(squeeze(signalArray(:, 2, :)), expected, ...
+ AbsTol=testCase.AbsTol)
+ end
+
+ function testImageMaskLeavesPixelsOutOfNeuropilSignal(testCase)
+ rois = createRois(testCase.ImageSize);
+ imageMask = true(testCase.ImageSize);
+ imageMask(:, 1:8) = false; % Columns to the left of roi 1
+ backgroundValue = 1000;
+ imageData = createTwoLevelStack(~imageMask, ...
+ testCase.NumFrames, backgroundValue, 9000);
+
+ signalArray = nansen.twophoton.roisignals.extractF( ...
+ imageData, rois, imageMask=imageMask, signalDataType="double");
+
+ expected = backgroundValue * ones(testCase.NumFrames, 3);
+ testCase.verifyEqual(squeeze(signalArray(:, 2, :)), expected, ...
+ AbsTol=testCase.AbsTol)
+ end
+
+ function testBatchExtractEqualsSerialExtract(testCase)
+ rois = createRois(testCase.ImageSize);
+ imageData = createRandomStack(testCase.ImageSize, testCase.NumFrames);
+
+ serialSignals = nansen.twophoton.roisignals.extractF(imageData, rois, ...
+ extractFcn=@nansen.twophoton.roisignals.extract.serialExtract, ...
+ roiMaskFormat="struct", signalDataType="double");
+ batchSignals = nansen.twophoton.roisignals.extractF(imageData, rois, ...
+ extractFcn=@nansen.twophoton.roisignals.extract.batchExtract, ...
+ roiMaskFormat="sparse", signalDataType="double");
+
+ testCase.verifyEqual(batchSignals, serialSignals, AbsTol=testCase.AbsTol)
+ end
+
+ function testSubsetOfRoisEqualsColumnsOfFullExtraction(testCase)
+ rois = createRois(testCase.ImageSize);
+ imageData = createRandomStack(testCase.ImageSize, testCase.NumFrames);
+
+ allSignals = nansen.twophoton.roisignals.extractF(imageData, rois);
+ selectedSignals = nansen.twophoton.roisignals.extractF( ...
+ imageData, rois, roiInd=[3, 1]);
+
+ testCase.verifyEqual(selectedSignals, allSignals(:, :, [3, 1]))
+ end
+
+ function testSecondOutputHoldsOptionsGivenByName(testCase)
+ rois = createRois(testCase.ImageSize);
+ imageData = createRandomStack(testCase.ImageSize, testCase.NumFrames);
+
+ [~, params] = nansen.twophoton.roisignals.extractF( ...
+ imageData, rois, createNeuropilMask=false, roiInd=[3, 1]);
+
+ testCase.verifyFalse(params.createNeuropilMask)
+ testCase.verifyEqual(params.roiInd, [3, 1])
+ end
+
+ function testSecondOutputHoldsResolvedDefaults(testCase)
+ rois = createRois(testCase.ImageSize);
+ imageData = createRandomStack(testCase.ImageSize, testCase.NumFrames);
+
+ [~, params] = nansen.twophoton.roisignals.extractF(imageData, rois);
+
+ testCase.verifyEqual(params.roiInd, 1:3)
+ testCase.verifyEqual(params.imageMask, true(testCase.ImageSize))
+ testCase.verifyEqual(params.extractFcn, ...
+ @nansen.twophoton.roisignals.extract.serialExtract)
+ testCase.verifyTrue(params.createNeuropilMask)
+ end
+
+ function testOptionsStructEqualsNameValueArguments(testCase)
+ rois = createRois(testCase.ImageSize);
+ imageData = createRandomStack(testCase.ImageSize, testCase.NumFrames);
+ optionsStruct = struct('createNeuropilMask', false, 'roiInd', 2);
+
+ [signalsFromStruct, paramsFromStruct] = ...
+ nansen.twophoton.roisignals.extractF(imageData, rois, optionsStruct);
+ [signalsFromNames, paramsFromNames] = ...
+ nansen.twophoton.roisignals.extractF(imageData, rois, ...
+ createNeuropilMask=false, roiInd=2);
+
+ testCase.verifySize(signalsFromStruct, [testCase.NumFrames, 1])
+ testCase.verifyEqual(signalsFromStruct, signalsFromNames)
+ testCase.verifyEqual(paramsFromStruct, paramsFromNames)
+ end
+
+ function testNameValueArgumentReplacesStructField(testCase)
+ rois = createRois(testCase.ImageSize);
+ imageData = createRandomStack(testCase.ImageSize, testCase.NumFrames);
+ optionsStruct = struct('createNeuropilMask', true);
+
+ signalArray = nansen.twophoton.roisignals.extractF( ...
+ imageData, rois, optionsStruct, createNeuropilMask=false);
+
+ testCase.verifySize(signalArray, [testCase.NumFrames, 1, 3])
+ end
+
+ function testDefaultParameterStructIsAccepted(testCase)
+ % The struct holds fields such as roiInd_ that configure the
+ % options editor and are not options of extractF.
+ rois = createRois(testCase.ImageSize);
+ imageData = createRandomStack(testCase.ImageSize, testCase.NumFrames);
+ optionsStruct = nansen.twophoton.roisignals.extract.getDefaultParameters();
+
+ signalsFromStruct = nansen.twophoton.roisignals.extractF( ...
+ imageData, rois, optionsStruct);
+ signalsFromDefaults = nansen.twophoton.roisignals.extractF(imageData, rois);
+
+ testCase.verifyEqual(signalsFromStruct, signalsFromDefaults)
+ end
+
+ function testInvalidValueInOptionsStructThrowsError(testCase)
+ rois = createRois(testCase.ImageSize);
+ imageData = createRandomStack(testCase.ImageSize, testCase.NumFrames);
+ optionsStruct = struct('pixelComputationMethod', 'max');
+
+ testCase.verifyError( ...
+ @() nansen.twophoton.roisignals.extractF(imageData, rois, optionsStruct), ...
+ "MATLAB:validators:mustBeMember")
+ end
+ end
+end
+
+function [rois, roiMasks] = createRois(imageSize)
+%createRois Two overlapping rois and one roi without neighbors
+ roiMasks = false([imageSize, 3]);
+ roiMasks(10:18, 10:18, 1) = true;
+ roiMasks(10:18, 16:24, 2) = true;
+ roiMasks(30:38, 40:48, 3) = true;
+
+ rois = [RoI('Mask', roiMasks(:, :, 1)), ...
+ RoI('Mask', roiMasks(:, :, 2)), ...
+ RoI('Mask', roiMasks(:, :, 3))];
+end
+
+function imageData = createRandomStack(imageSize, numFrames)
+%createRandomStack Stack of uniformly distributed uint16 pixel values
+ imageData = uint16(randi([100, 4000], [imageSize, numFrames]));
+end
+
+function imageData = createTwoLevelStack(isBright, numFrames, backgroundValue, brightValue)
+%createTwoLevelStack Stack with one value in bright pixels and one elsewhere
+ frame = backgroundValue * ones(size(isBright), 'single');
+ frame(isBright) = brightValue;
+ imageData = repmat(frame, 1, 1, numFrames);
+end
+
+function signals = meanInMasks(imageData, masks)
+%meanInMasks Mean pixel value per frame (rows) and mask (columns)
+ numFrames = size(imageData, 3);
+ numMasks = size(masks, 3);
+ pixelValues = reshape(double(imageData), [], numFrames);
+
+ signals = zeros(numFrames, numMasks);
+ for i = 1:numMasks
+ mask = masks(:, :, i);
+ signals(:, i) = mean(pixelValues(mask(:), :), 1)';
+ end
+end