diff --git a/.gitattributes b/.gitattributes new file mode 100644 index 0000000..5b68b4e --- /dev/null +++ b/.gitattributes @@ -0,0 +1,28 @@ +*.fig binary +*.mat binary +*.mdl binary diff merge=mlAutoMerge +*.mdlp binary +*.mexa64 binary +*.mexw64 binary +*.mexmaci64 binary +*.mlapp binary +*.mldatx binary +*.mlproj binary +*.mlx binary +*.p binary +*.sfx binary +*.sldd binary +*.slreqx binary merge=mlAutoMerge +*.slmx binary merge=mlAutoMerge +*.sltx binary +*.slxc binary +*.slx binary merge=mlAutoMerge +*.slxp binary + +## Other common binary file types +*.docx binary +*.exe binary +*.jpg binary +*.pdf binary +*.png binary +*.xlsx binary \ No newline at end of file diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..804e30f --- /dev/null +++ b/.gitignore @@ -0,0 +1,31 @@ +# Windows default autosave extension +*.asv + +# OSX / *nix default autosave extension +*.m~ + +# Compiled MEX binaries (all platforms) +*.mex* + +# Packaged app and toolbox files +*.mlappinstall +*.mltbx + +# Generated helpsearch folders +helpsearch*/ + +# Simulink code generation folders +slprj/ +sccprj/ + +# Matlab code generation folders +codegen/ + +# Simulink autosave extension +*.autosave + +# Simulink cache files +*.slxc + +# Octave session info +octave-workspace \ No newline at end of file diff --git a/datastore/DeepInterpolationDataStore.m b/datastore/DeepInterpolationDataStore.m new file mode 100644 index 0000000..aed9c4d --- /dev/null +++ b/datastore/DeepInterpolationDataStore.m @@ -0,0 +1,181 @@ +classdef DeepInterpolationDataStore < matlab.io.Datastore & ... + matlab.io.datastore.Subsettable & ... + matlab.io.datastore.Shuffleable + %DEEPINTERPOLATIONDATASTORE - custom datastore for DeepInterpolation-MATLAB + % + % This custom datastore provides flanking frames and the center frame + % upon one read. + % + % dids = DeepInterpolationDataStore(pathToTiffFile, options); + % + % Input must be full path to a grayscale, single channel tiff stack of + % >'flankingFrames + 2' frames in order to work (this is according to + % the function of the DeepInterpolation network). + % + % If the image frame dimensions differ from the network input size + % the doAutoResize option controls whether the frames will be + % automatically resized to outputFrameSize before they are returned. + % + % Dr. Thomas Künzel, The MathWorks, 2023 + % tkuenzel@mathworks.com + + properties + fileName string + doAutoResize logical + numberOfFlankingFrames int32 + outputFrameSize int32 + frameCount double + dsSetCount double + setsAvailable double + currentSet double + setFramesStartIndex double + currentFrameIndex double + allSetFrameStartIndices double + end + + methods % begin methods section + function myds = DeepInterpolationDataStore(filePath, options) + arguments + filePath + options = struct(); + end + if isfield(options, 'doAutoResize') + myds.doAutoResize = logical(options.doAutoResize); + else % default + myds.doAutoResize = false; + end + if isfield(options, 'numberOfFlankingFrames') + assert(options.numberOfFlankingFrames > 0) + assert(rem(options.numberOfFlankingFrames, 2) == 0, "Number of flanking frames should be even.") + myds.numberOfFlankingFrames = options.numberOfFlankingFrames; + else % default + myds.numberOfFlankingFrames = 60; + end + if isfield(options, 'outputFrameSize') + if not (myds.doAutoResize) + warning('Specified an output frame size, but auto-resize is turned off.') + end + outputFrameSize = options.outputFrameSize; + if not (isnumeric(outputFrameSize) && isequal(size(outputFrameSize), [1, 2]) && all(outputFrameSize >= 0)) + error('outputFrameSize should be an array with dimension 1x2 and non-negative values.'); + end + myds.outputFrameSize = outputFrameSize; + else % default + myds.outputFrameSize = [512, 512]; + end + assert(isfile(filePath)); + myds.fileName = filePath; + fileInfo = imfinfo(myds.fileName); + assert(strcmp(fileInfo(1).Format,'tif'), "Must be tiff format"); + if ~myds.doAutoResize + assert((fileInfo(1).Width == myds.outputFrameSize(1)) ... + && (fileInfo(1).Height == myds.outputFrameSize(2)),... + "Actual frame size is not equal to specified outputFrameSize"); + end + framePerSetCount = myds.numberOfFlankingFrames + 3; + assert(numel(fileInfo) >= framePerSetCount, "Not enough frames in stack for DeepInterpolation"); + myds.frameCount = numel(fileInfo); + myds.dsSetCount = myds.frameCount - framePerSetCount + 1; + myds.allSetFrameStartIndices = 1:myds.dsSetCount; + reset(myds); + end + + function tf = hasdata(myds) + % Return true if more data is available. + tf = myds.setsAvailable > 0; + end + + function [data,info] = read(myds) + % Read data and information about the extracted data. + assert(hasdata(myds), "No more data to read"); + rawRefFrame = imread(myds.fileName,myds.currentFrameIndex); + if myds.doAutoResize + rawRefFrame = imresize(rawRefFrame,myds.outputFrameSize); + end + refFrame = single(rawRefFrame); + flankingFrames = single( ... + ones(myds.outputFrameSize(1), ... + myds.outputFrameSize(2), ... + myds.numberOfFlankingFrames) ... + ); + framecount = 1; + for leftFrame = 0:(myds.numberOfFlankingFrames/2) - 1 + rawThisFrame = imread(myds.fileName,myds.setFramesStartIndex+leftFrame); + if myds.doAutoResize + rawThisFrame = imresize(rawThisFrame,myds.outputFrameSize); + end + thisFrame = single(rawThisFrame); + flankingFrames(:,:,framecount) = thisFrame; + framecount = framecount + 1; + end + %frame myds.numberOfFlankingFrames + 1 is center frame. One + %frame before and one after is left out. + startRightFrame = myds.numberOfFlankingFrames/2 + 3; + for rightFrame = startRightFrame:startRightFrame + myds.numberOfFlankingFrames/2 - 1 + rawThisFrame = imread(myds.fileName,myds.setFramesStartIndex+rightFrame); + if myds.doAutoResize + rawThisFrame = imresize(rawThisFrame,myds.outputFrameSize); + end + thisFrame = single(rawThisFrame); + flankingFrames(:,:,framecount) = thisFrame; + framecount = framecount + 1; + end + data = {flankingFrames, refFrame}; + info = "set=" + num2str(myds.currentSet) + "/centerframe=" + num2str(myds.currentFrameIndex); + myds.currentSet = myds.currentSet + 1; + myds.setsAvailable = myds.setsAvailable - 1; + if myds.setsAvailable > 0 + myds.setFramesStartIndex = myds.allSetFrameStartIndices(myds.currentSet); + myds.currentFrameIndex = myds.setFramesStartIndex + myds.numberOfFlankingFrames/2 + 1; + else + myds.setFramesStartIndex = nan; + myds.currentFrameIndex = nan; + end + end + + function reset(myds) + % Reset to the start of the data. + myds.setsAvailable = myds.dsSetCount; + myds.currentSet = 1; + myds.setFramesStartIndex = myds.allSetFrameStartIndices(myds.currentSet); + myds.currentFrameIndex = myds.setFramesStartIndex + myds.numberOfFlankingFrames/2 + 1; + end% + + function shmyds = shuffle(myds) + % shmyds = shuffle(myds) shuffles the sets + % Create a copy of datastore + shmyds = copy(myds); + + % Shuffle files and corresponding labels + numObservations = shmyds.dsSetCount; + idx = randperm(numObservations); + shmyds.allSetFrameStartIndices = myds.allSetFrameStartIndices(idx); + reset(shmyds); + end + end + + methods (Hidden = true) + function frac = progress(myds) + % Determine percentage of data read from datastore + if hasdata(myds) + frac = (myds.currentSet-1) / myds.dsSetCount; + else + frac = 1; + end + end + end + + methods (Access = protected) + function subds = subsetByReadIndices(myds, indices) + subds = copy(myds); + subds.dsSetCount = numel(indices); + subds.allSetFrameStartIndices = myds.allSetFrameStartIndices(indices); + reset(subds); + end + + function n = maxpartitions(myds) + n = myds.dsSetCount; + end + end + +end \ No newline at end of file diff --git a/datastore/test/TestDeepInterpolationDataStore.m b/datastore/test/TestDeepInterpolationDataStore.m new file mode 100644 index 0000000..da8d39f --- /dev/null +++ b/datastore/test/TestDeepInterpolationDataStore.m @@ -0,0 +1,390 @@ +% Testsuite for the DeepInterpolation Datastore +% https://www.mathworks.com/help/matlab/import_export/testing-guidelines-for-custom-datastores.html +% +% The test-suite expects the tiff-files "ophys_tiny_761605196.tif" and +% "ophys_tiny_761605196_wrongsize.tif" in the subfolder data of the current folder. +% If that is different, alter TESTSTACKFULLFILE! +% +% Note that a different test-dataset will fail many tests as they check for +% specific frames or framecounts. +% +% Dr. Thomas Künzel, The MathWorks, 2023 +% tkuenzel@mathworks.com + +base_folder = setup; + +TESTSTACKFULLFILE = fullfile(base_folder,"sample_data","ophys_tiny_761605196.tif"); +WRONGSIZETESTSTACKFULLFILE = fullfile(base_folder,"sample_data","ophys_tiny_761605196_wrongsize.tif"); + +%% Test 01: Check if matlab.io.Datastore is superclass. +% A deepinterpolation datastore inherits from abstract Datastore +dids = DeepInterpolationDataStore(TESTSTACKFULLFILE); +assert(isa(dids,'matlab.io.Datastore')); + +%% Test 02: Call the read +% A deepinterpolation datastore returns training data (512x512x60) and the reference +% frame (512x512) in a 1x2 cell array +% The first center- or reference frame to be read is actually the 32nd frame of the +% dataset, as deepinterpolation takes 60 flanking frames [n-31:n-2 n+2:n+31] around a +% center frame, leaving one frame each left and right as a gap, for training. +% The first frame that allows that is the 32nd! +% The combination of training data and center frame is called a "set" here, and the +% number of sets, following the logic above, is N-Frames - 62. +dids = DeepInterpolationDataStore(TESTSTACKFULLFILE); +t = read(dids); +assert(iscell(t)); +assert(numel(t) == 2); +assert(all((size(t{1})==[512 512 60]))); +assert(all((size(t{2})==[512 512]))); + +%% Test 03: Call the read again +% Standard sequential datastore behavior: advance to "second" center-frame +dids = DeepInterpolationDataStore(TESTSTACKFULLFILE); +read(dids); +[t,info] = read(dids); +assert(iscell(t)); +assert(numel(t) == 2); +assert(all((size(t{1})==[512 512 60]))); +assert(all((size(t{2})==[512 512]))); +assert(info == "set=2/centerframe=33") + +%% Test 04: Successfully read until no more data +% The testdatastack should return 38 sets of training frames and reference frame +dids = DeepInterpolationDataStore(TESTSTACKFULLFILE); +count = 0; +while(hasdata(dids)) + read(dids); + count = count + 1; +end +assert(count == 38); + +%% Test 05: Call read when out of data +dids = DeepInterpolationDataStore(TESTSTACKFULLFILE); +try + read(dids); +catch ME + assert(ME.message == "No more data to read"); +end + +%% Test 11: Call the readall method +% Standard sequential datastore behavior: return all data in a nx2 Cell array, for +% the test dataset n==38 +dids = DeepInterpolationDataStore(TESTSTACKFULLFILE); +allt = readall(dids); +assert(all(size(allt)==[38 2])); + +%% Test 12: Call the readall method when out of data +dids = DeepInterpolationDataStore(TESTSTACKFULLFILE); +while(hasdata(dids)) + read(dids); +end +allt = readall(dids); +assert(all(size(allt)==[38 2])); + +%% Test 21: Call the reset method before any read +% Reset to "first" frame (in fact, the 32nd frame in the stack...) +dids = DeepInterpolationDataStore(TESTSTACKFULLFILE); +reset(dids); +[t,info] = read(dids); +assert(iscell(t)); +assert(numel(t) == 2); +assert(all((size(t{1})==[512 512 60]))); +assert(all((size(t{2})==[512 512]))); +assert(info == "set=1/centerframe=32"); + +%% Test 22: Call the reset method after some reads +dids = DeepInterpolationDataStore(TESTSTACKFULLFILE); +read(dids); +read(dids); +read(dids); +reset(dids); +[t,info] = read(dids); +assert(iscell(t)); +assert(numel(t) == 2); +assert(all((size(t{1})==[512 512 60]))); +assert(all((size(t{2})==[512 512]))); +assert(info == "set=1/centerframe=32"); + +%% Test 23: Call the reset method after readall +dids = DeepInterpolationDataStore(TESTSTACKFULLFILE); +read(dids); +tall = readall(dids); +reset(dids); +[t,info] = read(dids); +assert(iscell(t)); +assert(numel(t) == 2); +assert(all((size(t{1})==[512 512 60]))); +assert(all((size(t{2})==[512 512]))); +assert(info == "set=1/centerframe=32"); + +%% Test 24: Call the reset method when out of data +dids = DeepInterpolationDataStore(TESTSTACKFULLFILE); +while(hasdata(dids)) + read(dids); +end +reset(dids); +[t,info] = read(dids); +assert(iscell(t)); +assert(numel(t) == 2); +assert(all((size(t{1})==[512 512 60]))); +assert(all((size(t{2})==[512 512]))); +assert(info == "set=1/centerframe=32"); + +%% Test 31: Call the progress method before any reads +dids = DeepInterpolationDataStore(TESTSTACKFULLFILE); +p = progress(dids); +assert(p==0); + +%% Test 32: Call the progress method after readall but before read. +dids = DeepInterpolationDataStore(TESTSTACKFULLFILE); +readall(dids); +p = progress(dids); +assert(p==0); + +%% Test 33: Call the progress method after read +dids = DeepInterpolationDataStore(TESTSTACKFULLFILE); +read(dids); +p = progress(dids); +assert((p-0.0270)<1e-4); + +%% Test 34: Call the progress method when out of data +dids = DeepInterpolationDataStore(TESTSTACKFULLFILE); +while(hasdata(dids)) + read(dids); +end +p = progress(dids); +assert((p-1.0)<1e-4); + +%% Test 35: Progress produces correct output over all sets +dids = DeepInterpolationDataStore(TESTSTACKFULLFILE); +idx = 1; +while(hasdata(dids)) + p(idx) = progress(dids); + read(dids); + idx = idx + 1; +end +p(idx) = progress(dids); +assert(all(p <= 1), "progress must be between 0 and 1!"); +assert(all(diff(p) > 0), "progress must be monotonically increasing"); + +%% Test 41: Preview methods returns first set +dids = DeepInterpolationDataStore(TESTSTACKFULLFILE); +t = preview(dids); +assert(iscell(t)); +assert(numel(t) == 2); +assert(all((size(t{1})==[512 512 60]))); +assert(all((size(t{2})==[512 512]))); +assert((t{1}(323,23,44)-138)<1e-4); + +%% Test 42: Preview methods returns first set after reads +dids = DeepInterpolationDataStore(TESTSTACKFULLFILE); +read(dids); +read(dids); +t = preview(dids); +assert(iscell(t)); +assert(numel(t) == 2); +assert(all((size(t{1})==[512 512 60]))); +assert(all((size(t{2})==[512 512]))); +assert((t{1}(323,23,44)-138)<1e-4); + +%% Test 43: Preview methods returns first set after readall +dids = DeepInterpolationDataStore(TESTSTACKFULLFILE); +readall(dids); +t = preview(dids); +assert(iscell(t)); +assert(numel(t) == 2); +assert(all((size(t{1})==[512 512 60]))); +assert(all((size(t{2})==[512 512]))); +assert((t{1}(323,23,44)-138)<1e-4); + +%% Test 44: Preview methods returns first set after reads and reset +dids = DeepInterpolationDataStore(TESTSTACKFULLFILE); +read(dids); +read(dids); +reset(dids); +t = preview(dids); +assert(iscell(t)); +assert(numel(t) == 2); +assert(all((size(t{1})==[512 512 60]))); +assert(all((size(t{2})==[512 512]))); +assert((t{1}(323,23,44)-138)<1e-4); + +%% Test 45: Preview methods returns first set when out of data +dids = DeepInterpolationDataStore(TESTSTACKFULLFILE); +while hasdata(dids) + read(dids); +end +t = preview(dids); +assert(iscell(t)); +assert(numel(t) == 2); +assert(all((size(t{1})==[512 512 60]))); +assert(all((size(t{2})==[512 512]))); +assert((t{1}(323,23,44)-138)<1e-4); + +%% Test 46: Read after preview is from correct location +dids = DeepInterpolationDataStore(TESTSTACKFULLFILE); +read(dids); +read(dids); +preview(dids); +[t,info] = read(dids); +assert(iscell(t)); +assert(numel(t) == 2); +assert(all((size(t{1})==[512 512 60]))); +assert(all((size(t{2})==[512 512]))); +assert(info == "set=3/centerframe=34"); + +%% Test 47: Readall after preview works correctly +dids = DeepInterpolationDataStore(TESTSTACKFULLFILE); +preview(dids); +allt = readall(dids); +assert(all(size(allt)==[38 2])); + +%% Test 48: Hasdata works correctly after preview +dids = DeepInterpolationDataStore(TESTSTACKFULLFILE); +read(dids); +preview(dids); +assert(hasdata(dids)); + +%% Test 50: Call Partition, check if the partition is a deepinterpolation datastore +% A partition of the deepinterpolation datastore is merely pointing to a different +% startframe and has a smaller number of available sets. All the other behaviors are +% derived from that +dids = DeepInterpolationDataStore(TESTSTACKFULLFILE); +subds = partition(dids,10,7); +assert(isa(subds,'matlab.io.Datastore')); + +%% Test 51: Call Partition, check if the partition has correct properties +dids = DeepInterpolationDataStore(TESTSTACKFULLFILE); +subds = partition(dids,10,7); +assert(isequal(properties(dids),properties(subds))); + +%% Test 52: Call Partition, check if the partition has correct methods +dids = DeepInterpolationDataStore(TESTSTACKFULLFILE); +subds = partition(dids,10,7); +assert(isequal(methods(dids),methods(subds))); + +%% Test 53: Call Partition, successfully read from a partition +dids = DeepInterpolationDataStore(TESTSTACKFULLFILE); +subds = partition(dids,10,7); +[t,info] = read(subds); +assert(iscell(t)); +assert(numel(t) == 2); +assert(all((size(t{1})==[512 512 60]))); +assert(all((size(t{2})==[512 512]))); +assert(info == "set=1/centerframe=56"); + +%% Test 54: Call last Partition, read all data +dids = DeepInterpolationDataStore(TESTSTACKFULLFILE); +subds = partition(dids,13,13); +while hasdata(subds) + read(subds); +end + +%% Test 55: Call Partition, check that total number of sets/frames is correct +dids = DeepInterpolationDataStore(TESTSTACKFULLFILE); +partsets = 0; +for iii = 1:10 + subds = partition(dids,10,iii); + partsets = partsets + subds.dsSetCount; +end +assert(partsets-38 < 1e-4); + +%% Test 56: Call partition with partition=1 and index=1 and validate +dids = DeepInterpolationDataStore(TESTSTACKFULLFILE); +subds = partition(dids,1,1); +assert(isequal(properties(dids),properties(subds))); +assert(isequal(methods(dids),methods(subds))); +assert(isequaln(read(subds),read(dids))); +assert(isequaln(preview(subds),preview(dids))); + +%% Test 57: Partition a partition again and validate +dids = DeepInterpolationDataStore(TESTSTACKFULLFILE); +subds = partition(dids,2,2); +subsubds = partition(subds,2,2); +[t,info] = read(subsubds); +assert(iscell(t)); +assert(numel(t) == 2); +assert(all((size(t{1})==[512 512 60]))); +assert(all((size(t{2})==[512 512]))); +assert(info == "set=1/centerframe=61"); + +%% Test 58: Progress is working correctly in partitions +dids = DeepInterpolationDataStore(TESTSTACKFULLFILE); +subds = partition(dids,2,2); +idx = 1; +clear p +while(hasdata(subds)) + p(idx) = progress(subds); + read(subds); + idx = idx + 1; +end +p(idx) = progress(subds); +assert(p(1)<1e-4, "first progress item should be 0"); +assert((p(end)-1)<1e-5, "last progress item should be 1"); +assert(all(p <= 1), "progress must be between 0 and 1!"); +assert(all(diff(p) > 0), "progress must be monotonically increasing"); + +%% Test 60: Create a Subset of the datastore and validate +% deepinterpolation datastore is subsettable, again this just assigns the correct +% startframe and number of available sets, but points to the same file +dids = DeepInterpolationDataStore(TESTSTACKFULLFILE); +subds = subset(dids,5:20); +[t,info]=read(subds); +assert(iscell(t)); +assert(numel(t) == 2); +assert(all((size(t{1})==[512 512 60]))); +assert(all((size(t{2})==[512 512]))); +assert(info == "set=1/centerframe=36"); + +%% Test 61: Progress is working correctly in subsets +dids = DeepInterpolationDataStore(TESTSTACKFULLFILE); +subds = subset(dids,5:10); +idx = 1; +clear p +while(hasdata(subds)) + p(idx) = progress(subds); + read(subds); + idx = idx + 1; +end +p(idx) = progress(subds); +assert(p(1)<1e-4, "first progress item should be 0"); +assert((p(end)-1)<1e-5, "last progress item should be 1"); +assert(all(p <= 1), "progress must be between 0 and 1!"); +assert(all(diff(p) > 0), "progress must be monotonically increasing"); + +%% Test 71: Fail on wrong image size +try + dids = DeepInterpolationDataStore(WRONGSIZETESTSTACKFULLFILE); %#ok +catch EM + assert(strcmp(EM.message,'Actual frame size is not equal to specified outputFrameSize'),"Fail on wrong tiff size"); +end + +%% Test 72: Read data with automatic resize +options.doAutoResize = true; +dids = DeepInterpolationDataStore(WRONGSIZETESTSTACKFULLFILE, options); +t = read(dids); +assert(iscell(t)); +assert(numel(t) == 2); +assert(all((size(t{1})==[512 512 60]))); +assert(all((size(t{2})==[512 512]))); + +%% Test 81: Shuffling returns randomized datastore without error +dids = DeepInterpolationDataStore(TESTSTACKFULLFILE); +shds = shuffle(dids); +assert(isa(shds,'matlab.io.Datastore')); + +%% Test 82: can read sequentially from shuffled datastore +dids = DeepInterpolationDataStore(TESTSTACKFULLFILE); +shds = shuffle(dids); +count = 0; +while(hasdata(shds)) + read(shds); + count = count + 1; +end +assert(count == 38); + +%% Test 83: can read all from shuffled datastore +dids = DeepInterpolationDataStore(TESTSTACKFULLFILE); +shds = shuffle(dids); +allt = readall(shds); +assert(all(size(allt)==[38 2])); \ No newline at end of file diff --git a/examples/customdatastore_example.mlx b/examples/customdatastore_example.mlx new file mode 100644 index 0000000..b6598e1 Binary files /dev/null and b/examples/customdatastore_example.mlx differ diff --git a/sample_data/denoised_ophys_tiny_761605196.tif b/sample_data/denoised_ophys_tiny_761605196.tif new file mode 100644 index 0000000..4663c07 Binary files /dev/null and b/sample_data/denoised_ophys_tiny_761605196.tif differ diff --git a/sample_data/ophys_tiny_761605196_wrongsize.tif b/sample_data/ophys_tiny_761605196_wrongsize.tif new file mode 100644 index 0000000..e1e45b4 Binary files /dev/null and b/sample_data/ophys_tiny_761605196_wrongsize.tif differ diff --git a/sample_data/pretrainedNetwork.mat b/sample_data/pretrainedNetwork.mat new file mode 100644 index 0000000..27ed6fc Binary files /dev/null and b/sample_data/pretrainedNetwork.mat differ diff --git a/sample_data/retrainedNetwork_20230305_195921.mat b/sample_data/retrainedNetwork_20230305_195921.mat new file mode 100644 index 0000000..d7d67e7 Binary files /dev/null and b/sample_data/retrainedNetwork_20230305_195921.mat differ diff --git a/setup.m b/setup.m index b83456e..ec57e35 100644 --- a/setup.m +++ b/setup.m @@ -5,4 +5,5 @@ addpath(fullfile(base_folder, "network_layers")) addpath(fullfile(base_folder, "sample_data")) addpath(fullfile(base_folder, "examples")) + addpath(fullfile(base_folder, "datastore")) end \ No newline at end of file