windows 10 test OK

README.md

+This code has been tested under Windows 10 version 1903 using Matlab 2019b with all toolboxes installed.
+To run this code example you need to keep the folder structure provided in this repository as depicted below:
+
+CloneClassifier
+|
+|_dataIn
+|
+|_dataOut
+|
+|_src
+
+The image data folder needs to be placed into the folder dataIn
+
+To complete the dependencies, please add the following files to the src folder:
+
+-dirrec.m (can be downloaded from 'https://nl.mathworks.com/matlabcentral/fileexchange/15505-recursive-dir')
+-xmlreadstring.m (can be downloaded from 'https://nl.mathworks.com/matlabcentral/fileexchange/46844-deployanywhere-zip')
+
+To trigger code execution run the main.m script located in src from Matlab
\ No newline at end of file

src/dirrec.m

+function [varargout] = dirrec(reper,ext)
+%
+% Find files recursively in a given folder.
+%
+%     C=dirrec('c:\windows') returns a cell C with the full pathname of all
+%      files in the c:\windows folder and all its sub-folders.
+%
+%     C=dirrec('c:\windows','.exe') idem but returns only the files with
+%      extension .exe.
+%
+%     C=dirrec('c:\windows','co*') idem but returns only the files starting with
+%      the two letters co (comsetup.log, control.ini, ...).
+%
+%     C=dirrec('c:\windows',{'.exe','.dll') idem but returns files with both
+%      .exe and .dll extensions.
+%
+%     dirrec('c:\windows','.cmd') only displays the list of the .cmd files in
+%      the Matlab command window
+%
+%           c:\windows\system32\login.cmd
+%           c:\windows\system32\usrlogon.cmd
+%
+%     Note that extension should be given in lower case.
+%
+%     See also DIR.
+%
+%     Luc Masset (2007)
+
+%  Algorithm:
+% The three cases - no extension, one extension and multiple extensions - are separeted to speed up the
+% search process. The function fileextDR replaces the fileparts function because we only need the
+% extension.
+
+%initialisation
+if nargout,
+ varargout=[];
+end
+listF=[];
+
+%input arguments
+if ~nargin | nargin > 2,
+ error('DIRREC requires 1 or two arguments')
+ return
+elseif nargin == 1,
+ ext=[];
+else
+ if ~iscell(ext),
+  if strcmpi(ext,'.*'),
+   ext=[];
+  end
+ end
+end
+
+%list of folders
+listD{1,1}=reper;       % a cell containing all the searched folders
+indD(1)=1;              % a vector (same size as listD) indicating
+                        % that a folder has been searched (1) or not (0)
+
+%cases
+if isempty(ext),
+ icase=1;
+elseif ~iscell(ext),
+ ii=strfind(ext,'*');
+ if length(ii) == 1 & ii == length(ext);
+  icase=2;
+ else
+  icase=3;
+ end
+elseif iscell(ext),
+ icase=4;
+else
+ ext=[];
+ icase=1;
+end
+
+%case 1: no extension given
+switch icase
+case 1,
+ while 1,
+  ind=find(indD);
+  if isempty(ind),
+   break;
+  end
+  ind=ind(1);
+  rep=listD{ind};
+  [listdir,listfile]=getdirDR1(rep);
+  listF=[listF listfile];
+  indD(ind)=0;
+  nbd=length(listdir);
+  if nbd,
+   listD=[listD listdir];
+   indD=[indD ones(1,nbd)];
+  end
+ end
+
+%case 2: start of the name given (--*)
+case 2,
+ stc=strrep(ext,'*','');
+ while 1,
+  ind=find(indD);
+  if isempty(ind),
+   break;
+  end
+  ind=ind(1);
+  rep=listD{ind};
+  [listdir,listfile]=getdirDR4(rep,stc);
+  listF=[listF listfile];
+  indD(ind)=0;
+  nbd=length(listdir);
+  if nbd,
+   listD=[listD listdir];
+   indD=[indD ones(1,nbd)];
+  end
+ end
+
+%case 3: only one extension given
+case 3,
+ while 1,
+  ind=find(indD);
+  if isempty(ind),
+   break;
+  end
+  ind=ind(1);
+  rep=listD{ind};
+  [listdir,listfile]=getdirDR2(rep,ext);
+  listF=[listF listfile];
+  indD(ind)=0;
+  nbd=length(listdir);
+  if nbd,
+   listD=[listD listdir];
+   indD=[indD ones(1,nbd)];
+  end
+ end
+
+%case 4: several extensions given
+case 4,
+ while 1,
+  ind=find(indD);
+  if isempty(ind),
+   break;
+  end
+  ind=ind(1);
+  rep=listD{ind};
+  [listdir,listfile]=getdirDR3(rep,ext);
+  listF=[listF listfile];
+  indD(ind)=0;
+  nbd=length(listdir);
+  if nbd,
+   listD=[listD listdir];
+   indD=[indD ones(1,nbd)];
+  end
+ end
+end
+
+%display results
+if ~nargout,
+ for i=1:length(listF),
+  fprintf('%s\n',listF{i});
+ end
+else
+ varargout{1}=listF;
+end
+
+return
+
+%------------------------------------------------------------------------------
+function [listdir,listfile] = getdirDR1(reper)
+
+%dir of folder reper
+S=dir(reper);
+
+%separate sub-folders of reper and files
+n=size(S,1);
+listdir=cell(1,n);      % list of sub-folders
+listfile=cell(1,n);     % list of files
+for i=1:n,
+ name=S(i).name;
+ if S(i).isdir,
+  if strcmp(name,'.'),  % remove current folder (.)
+   continue;
+  end
+  if strcmp(name,'..'), % remove parent folder (..)
+   continue;
+  end
+  listdir{i}=fullfile(reper,S(i).name);
+ else
+  listfile{i}=fullfile(reper,S(i).name);
+ end
+end
+
+%reorder results
+ind=find(cellfun('isempty',listdir));
+listdir(ind)=[];
+ind=find(cellfun('isempty',listfile));
+listfile(ind)=[];
+
+return
+
+%------------------------------------------------------------------------------
+function [listdir,listfile] = getdirDR2(reper,ext)
+
+%dir of folder reper
+S=dir(reper);
+
+%separate sub-folders of reper and files
+n=size(S,1);
+listdir=cell(1,n);      % list of sub-folders
+listfile=cell(1,n);     % list of files
+nd=0;
+nf=0;
+for i=1:n,
+ name=S(i).name;
+ if S(i).isdir,
+  if strcmp(name,'.'),  % remove current folder (.)
+   continue;
+  end
+  if strcmp(name,'..'), % remove parent folder (..)
+   continue;
+  end
+  nd=nd+1;
+  listdir{nd}=fullfile(reper,S(i).name);
+ else
+  exte=fileextDR(name);       % compare extension
+  if strcmpi(exte,ext),       % with given extension
+   nf=nf+1;
+   listfile{nf}=fullfile(reper,S(i).name);
+  end
+ end
+end
+
+%reorder results
+listdir(nd+1:end)=[];
+listfile(nf+1:end)=[];
+
+return
+
+%------------------------------------------------------------------------------
+function [listdir,listfile] = getdirDR3(reper,ext)
+
+%dir of folder reper
+S=dir(reper);
+
+%separate sub-folders of reper and files
+n=size(S,1);
+listdir=cell(1,n);      % list of sub-folders
+listfile=cell(1,n);     % list of files
+for i=1:n,
+ name=S(i).name;
+ if S(i).isdir,
+  if strcmp(name,'.'),  % remove current folder (.)
+   continue;
+  end
+  if strcmp(name,'..'), % remove parent folder (..)
+   continue;
+  end
+  listdir{i}=fullfile(reper,S(i).name);
+ else
+  exte=fileextDR(name);       % extension of the file
+  if isempty(exte),
+   continue;
+  end
+  if strmatch(lower(exte),ext,'exact'),   % compare extension with given extensions
+   listfile{i}=fullfile(reper,S(i).name);
+  end
+ end
+end
+
+%reorder results
+ind=find(cellfun('isempty',listdir));
+listdir(ind)=[];
+ind=find(cellfun('isempty',listfile));
+listfile(ind)=[];
+
+return
+
+%------------------------------------------------------------------------------
+function [listdir,listfile] = getdirDR4(reper,stc)
+
+%dir of folder reper
+S=dir(reper);
+
+%separate sub-folders of reper and files
+n=size(S,1);
+listdir=cell(1,n);      % list of sub-folders
+listfile=cell(1,n);     % list of files
+nd=0;
+nf=0;
+for i=1:n,
+ name=S(i).name;
+ if S(i).isdir,
+  if strcmp(name,'.'),  % remove current folder (.)
+   continue;
+  end
+  if strcmp(name,'..'), % remove parent folder (..)
+   continue;
+  end
+  nd=nd+1;
+  listdir{nd}=fullfile(reper,S(i).name);
+ else
+  ii=strfind(name,stc);
+  if isempty(ii) | ii ~= 1,
+   continue;
+  end
+  nf=nf+1;
+  listfile{nf}=fullfile(reper,S(i).name);
+ end
+end
+
+%reorder results
+listdir(nd+1:end)=[];
+listfile(nf+1:end)=[];
+
+return
+
+%------------------------------------------------------------------------------
+function [ext] = fileextDR(fname)
+
+ext=[];
+ind=strfind(fname,filesep);
+if ~isempty(ind),
+ fname=fname(max(ind)+1:end);
+end
+ind=strfind(fname,'.');
+if isempty(ind),
+ return
+end
+ext=fname(max(ind):end);
+
+return

src/f_CloneClassifier.m

+function [Objects, Preview] = f_CloneClassifier(row, col, ch1Green, ch2Red, ch3Blue, SavePath, ShrinkFactor, Barcode)
+%Classify clones and create previews for a plate overview
+%   Detailed explanation goes here
+
+    %% Image analysis
+
+    % Red pos
+    %RedLP = imfilter(ch1Green, fspecial('gaussian', 60, 20), 'symmetric');
+    RedLP = imfilter(ch2Red, fspecial('gaussian', 21, 7), 'symmetric');
+    %it(ch2Red)
+    %it(RedLP)
+    RedPosMask = RedLP > 125;
+    %it(RedPosMask)
+
+    % Green Pos
+    %GreenLP = imfilter(ch2Red, fspecial('gaussian', 21, 7), 'symmetric');
+    GreenLP = imfilter(ch1Green, fspecial('gaussian', 60, 20), 'symmetric');
+    %it(GreenLP)
+    %GreenPosMask = GreenLP > 150;
+    GreenPosMask = GreenLP > 200;
+    %it(GreenPosMask)
+
+    SamplingZone = RedPosMask | GreenPosMask;
+    %it(SamplingZone)
+    
+    BlueLP = imfilter(ch3Blue, fspecial('gaussian', 60, 20), 'symmetric');
+    BluePosMask = BlueLP > 200;
+
+    %% Features
+
+
+
+    GreenVec = ch2Red(SamplingZone);
+    RedVec = ch1Green(SamplingZone);
+    Objects = table();
+    Objects.Area = sum(SamplingZone(:));
+    Objects.AreaRedPos = sum(RedPosMask(:));
+    Objects.AreaGreenPos = sum(GreenPosMask(:));
+    Objects.AreaBluePos = sum(BluePosMask(:));
+    Objects.GreenMean = mean(GreenVec);
+    Objects.RedMean = mean(RedVec);
+    Objects.Barcode = Barcode;
+    Objects.ROW = row;
+    Objects.COL = col;
+
+%     if (Objects.AreaRedPos / Objects.Area > 0.5) & ~(Objects.AreaGreenPos / Objects.Area > 0.5) & (Objects.Area >= 1e5)
+%         Class = 'Red';
+%     elseif ~(Objects.AreaRedPos / Objects.Area > 0.5) & (Objects.AreaGreenPos / Objects.Area > 0.5) & (Objects.Area >= 1e5)
+%         Class = 'Green';
+%     elseif (Objects.AreaRedPos / Objects.Area > 0.5) & (Objects.AreaGreenPos / Objects.Area > 0.5) & (Objects.Area >= 1e5)
+%         Class = 'RedGreen';
+%     else
+%         Class = 'Neg';
+%     end
+    if Objects.AreaBluePos / Objects.Area > 0.5
+        Class = 'Blue';
+    else
+        if (Objects.AreaRedPos / Objects.Area > 0.9) & ~(Objects.AreaGreenPos / Objects.Area > 0.1) & (Objects.Area >= 2e4)
+            Class = 'Red';
+        elseif ~(Objects.AreaRedPos / Objects.Area > 0.1) & (Objects.AreaGreenPos / Objects.Area > 0.9) & (Objects.Area >= 2e4)
+            Class = 'Green';
+        elseif (Objects.AreaRedPos / Objects.Area > 0.1) & (Objects.AreaGreenPos / Objects.Area > 0.1) & (Objects.Area >= 2e4)
+            Class = 'RedGreen';
+        else
+            Class = 'Neg';
+        end
+    end
+
+    Objects.Class = {Class};
+
+    %% Previews
+
+    %RGB = cat(3, imadjust(ch2Red, [0 0.01], [0 1]), imadjust(ch1Green, [0 0.01], [0 1]), zeros(size(ch1Green), 'uint16'));
+    RGB = cat(3, imadjust(ch2Red, [0 0.01], [0 1]), imadjust(ch1Green, [0 0.01], [0 1]),  imadjust(ch3Blue, [0 0.02], [0 1]));
+    % imtool(RGB)
+    Preview = imoverlay2(RGB, imdilate(bwperim(SamplingZone), strel('disk', 10)), [1 1 1]);
+    Preview = insertText(Preview, [1 1], Class, 'FontSize', 200, 'BoxColor', 'black', 'TextColor', 'white');
+    Preview = imresize(Preview, 1/ShrinkFactor);
+    % imtool(Preview)
+
+end
+

src/f_Local_Read_Files_Timeseries_Mosaic_BT2.m

+function [files, info, areaNames, plateLayout, areaMap] = f_Local_Read_Files_Timeseries_Mosaic_BT2(hcsPath, barcode, timepoint)
+%Author: Paul Antony 20141024
+%Loads image data and metadata from a time series barcode from the hcs drive
+%   hcsPath: Full path to the folder containing the meas folder, e.g., 'S:\OperaQEHS\OperaArchiveCol'
+%   barcode: The barcode to be analyzed
+%   timepoint: the counter of the MEAS folder to be analyzed
+%   Example: [files, info, areaNames, plateLayout, areaMap] = f_Local_Read_Files_Timeseries('Y', 'OB 20141022_MDPD1_3Stains_2Controls', 1)
+
+TimepointsStruct = dir([hcsPath, filesep, barcode]);
+TimepointsTable = struct2table(TimepointsStruct);
+TimepointsTableClean = varfun(@(x) regexp(x,'Meas'), TimepointsTable, 'InputVariables', {'name'});
+TimepointsTableFilter = cellfun(@(x) ~isempty(x), TimepointsTableClean.Fun_name, 'UniformOutput', false);
+TimepointsTableFilter = cell2mat(TimepointsTableFilter);
+TimepointsTableClean2 = TimepointsTable(TimepointsTableFilter, TimepointsTable.Properties.VariableNames);
+TimepointFolders = TimepointsTableClean2.name;
+%files = dirrec([hcsDrive, ':\OperaArchiveCol\OB 20141022_MDPD1_3Stains_2Controls', '\', TimepointFolders{timepoint}], '.flex');
+files = dirrec([hcsPath, filesep, barcode, filesep, TimepointFolders{timepoint}], '.flex');
+%[info, areaNames, plateLayout, areaMap] = flexinfo(files);
+[info, areaNames, plateLayout, areaMap] = flexinfoMosaic(files);
+
+end
+
+
+
+

src/f_LogDependencies.m

+function [FileList] = f_LogDependencies(Name, SavePath)
+%Collect all upstream scripts and save them in a folder named Dependencies
+%   Name: Function or script name as string
+%   SavePath: 
+%   Example: f_LogDependencies('Workflow_RosellaLC3_BPF_Summary_20170313.m', 'S:\HCS_Platform\Data\JavierJarazo\Rosella\LC3\Analysis_20170323_143443')
+    
+    SavePath = [SavePath filesep 'Dependencies'];
+    mkdir(SavePath);
+    
+    %% Scan recursively %%
+    FileList = matlab.codetools.requiredFilesAndProducts(Name)';
+    %%%%%%%%%%%%%%%%%%%%%%
+    FileCount = size(FileList,1);
+
+    for i = 1:FileCount
+        FileThis = regexp(FileList{i}, '\\.*\\(.*)', 'tokens'); FileThis = FileThis{:}{:};
+        SavePathThisFile = [SavePath, filesep, FileThis];
+        copyfile(FileList{i}, SavePathThisFile);
+    end
+    
+end
+

src/f_load_mosaic_image3D_Sparse.m

+function [XYmosaic, XYmosaicCell, XPositions, YPositions, ResolutionXY] = f_load_mosaic_image3D_Sparse(files, InfoTable, row, column, channel)
+%UNTITLED2 Summary of this function goes here
+%   Detailed explanation goes here
+        if isfield(table2struct(InfoTable), 'row') % For compatibility with data loading before 20161201
+            InfoRowFilter = (InfoTable.row == row) & (InfoTable.column == column);
+        elseif isfield(table2struct(InfoTable), 'Row') % Recommended use
+            InfoRowFilter = (InfoTable.Row == row) & (InfoTable.Column == column);
+        end
+        InfoSelection = InfoTable(InfoRowFilter, InfoTable.Properties.VariableNames);
+        CoordinateStrings = [InfoSelection.PositionX, InfoSelection.PositionY];
+        CoordinatesNumeric = cellfun(@(x) str2double(x), CoordinateStrings);
+        InfoSelection.PositionX = CoordinatesNumeric(:,1);
+        InfoSelection.PositionY = CoordinatesNumeric(:,2);
+
+        InfoSelectionSorted = sortrows(InfoSelection, {'PositionY', 'PositionX'}, 'descend');
+        
+        XPositions = unique(InfoSelectionSorted.PositionX);
+        YPositions = unique(InfoSelectionSorted.PositionY);
+        ResolutionXY = unique(InfoSelectionSorted.ResolutionXY);
+        XYmosaicCell = {};
+        
+        for f = 1:max(InfoSelectionSorted.field(:)) % loop over fields
+            InfoThisField = InfoSelectionSorted((InfoSelectionSorted.field == f), :);
+            cube = readflexcube(InfoThisField.filename{:} , 'PlaneCount', 1);
+            ImageThisField = cube.data(:, :, 1, [channel:InfoThisField.ChannelCount:(InfoThisField.ChannelCount*InfoThisField.PlaneCount)]); %Hoechst
+            % it(ImageThisField);
+            rowThisField = find(YPositions == InfoThisField.PositionY);
+            columnThisField = find(XPositions == InfoThisField.PositionX);
+            XYmosaicCell{rowThisField, columnThisField} = ImageThisField;
+        end
+        
+        EmptyCells = cellfun(@(x) isempty(x), XYmosaicCell);
+        ImSize = size(XYmosaicCell{min(find(~EmptyCells))});
+        XYmosaicCell(EmptyCells) = {zeros(ImSize, 'uint16')};
+        XYmosaic = cell2mat(XYmosaicCell);
+        % it(XYmosaic)
+
+end
+
+

src/flexfileinfoMosaic.m

+function [flexinfo, info] = flexfileinfoMosaic(fname)
+
+flexinfo.filename = fname;
+
+retries = 5;
+while retries > 0, try %#ok<ALIGN>
+    info = imfinfo(fname);
+    retries = 0;
+catch err
+    retries = retries - 1;
+    if retries == 0
+        rethrow(err)
+    end
+    fprintf('Unable to open file %s, %d retries left...\n', fname, retries);
+end, end
+
+xmlstr = info(1).UnknownTags(1).Value;
+xml = xmlreadstring(xmlstr);
+flexinfo.areaname = char(xml.getElementsByTagName('AreaName').item(0).getTextContent());
+flexinfo.barcode = char(xml.getElementsByTagName('Barcode').item(0).getTextContent());
+attr = xml.getElementsByTagName('WellCoordinate').item(0).getAttributes();
+flexinfo.row = str2double(char(attr.getNamedItem('Row').getNodeValue()));
+flexinfo.column = str2double(char(attr.getNamedItem('Col').getNodeValue()));
+flexinfo.meas = str2double(char(xml.getElementsByTagName('OperaMeasNo').item(0).getTextContent()));
+flexinfo.field = str2double(char( xml.getElementsByTagName('Images').item(0).getElementsByTagName('Sublayout').item(0).getTextContent()));
+flexinfo.timestamp = char(xml.getElementsByTagName('DateTime').item(0).getTextContent());
+flexinfo.PositionX = char(xml.getElementsByTagName('PositionX').item(0).getTextContent());
+flexinfo.PositionY = char(xml.getElementsByTagName('PositionY').item(0).getTextContent());
+flexinfo.PositionZ = char(xml.getElementsByTagName('PositionZ').item(0).getTextContent());
+Channels = {};
+for i = 0:xml.getElementsByTagName('Image').getLength-1
+    Channels{i+1} = char(xml.getElementsByTagName('Image').item(i).getElementsByTagName('CameraRef').item(0).getTextContent());
+end
+flexinfo.ChannelCount = length(unique(Channels));
+flexinfo.PlaneCount = xml.getElementsByTagName('Image').getLength / flexinfo.ChannelCount;
+flexinfo.ResolutionXY = eval(char(xml.getElementsByTagName('Image').item(i).getElementsByTagName('ImageResolutionX').item(0).getTextContent()));
+
+
+end
+

src/flexinfoMosaic.m

+function [info, areanames, areas, map] = flexinfoMosaic(files)
+% FLEXINFO extract metadata from flex files
+% [info, areaNames, plateLayout, areaMap] = FLEXFILES(wildcard) 
+% extracts metadata from the flex files matching wildcard
+% flexinfo is a struct with the following fields: 
+% - filename: full path to the flex file
+% - row, column: position of the flex file in the plate
+% - area: numerical identifier for the plate area
+% - areaname: name of the area from the plate layout
+% areaNames is a cell array with areanames indexed by area ids
+% plateLayout is a matrix with the assignment of area ids to well positions
+% areaMap is a containers.Map mapping area names to area ids
+
+%[~, files] = fileattrib(wildcard);
+%info = struct([]);
+%parfor i=1:length(files)
+for i=1:length(files)
+    info(i) = flexfileinfoMosaic(files{i});
+end
+
+[rows, columns] = flexplateinfo(info(1).filename);
+areas = nan(rows, columns);
+
+map = containers.Map();
+
+k = 1;
+
+areanames = cell(1);
+
+for i=1:length(info)
+    row = info(i).row;
+    column = info(i).column;
+    areaname = info(i).areaname;
+    if isKey(map, areaname)
+        area = map(areaname);
+        areas(row, column) = area;
+        info(i).area = area;
+    else
+        areas(row, column) = k;
+        map(areaname) = k;
+        info(i).area = k;
+        areanames{k} = areaname;
+        k = k+1;
+    end
+end
+
+end
+

src/flexplateinfo.m

+function [rows, cols] = flexplateinfo(fname)
+
+info = imfinfo(fname);
+xmlstr = info(1).UnknownTags(1).Value;
+xml = xmlreadstring(xmlstr);
+cols = str2double(char(xml.getElementsByTagName('YSize').item(0).getTextContent()));
+rows = str2double(char(xml.getElementsByTagName('XSize').item(0).getTextContent()));
+
+end
+

src/imoverlay2.m

+function out = imoverlay2(in, mask, color)
+%IMOVERLAY Create a mask-based image overlay.
+%   OUT = IMOVERLAY(IN, MASK, COLOR) takes an input image, IN, and a binary
+%   image, MASK, and produces an output image whose pixels in the MASK
+%   locations have the specified COLOR.
+%
+%   IN should be a grayscale or an RGB image of class uint8, uint16, int16,
+%   logical, double, or single.  If IN is double or single, it should be in
+%   the range [0, 1].  If it is not in that range, you might want to use
+%   mat2gray to scale it into that range.
+%
+%   MASK should be a two-dimensional logical matrix.
+%
+%   COLOR should be a 1-by-3 vector of values in the range [0, 1].  [0 0 0]
+%   is black, and [1 1 1] is white.
+%
+%   OUT is a uint8 RGB image.
+%
+%   Examples
+%   --------
+%   Overlay edge detection result in green over the original image.
+%       
+%       I = imread('cameraman.tif');
+%       bw = edge(I, 'canny');
+%       rgb = imoverlay(I, bw, [0 1 0]);
+%       imshow(rgb)
+%
+%   Treating the output of peaks as an image, overlay the values greater than
+%   7 in red.  The output of peaks is not in the usual grayscale image range
+%   of [0, 1], so use mat2gray to scale it.
+%
+%       I = peaks;
+%       mask = I > 7;
+%       rgb = imoverlay(mat2gray(I), mask, [1 0 0]);
+%       imshow(rgb, 'InitialMagnification', 'fit')
+
+%   Steven L. Eddins
+%   Copyright 2006-2012 The MathWorks, Inc.
+
+% If the user doesn't specify the color, use white.
+DEFAULT_COLOR = [1 1 1];
+if nargin < 3
+    color = DEFAULT_COLOR;
+end
+
+% Force the 2nd input to be logical.
+mask = (mask ~= 0);
+
+% Make the uint8 the working data class.  The output is also uint8.
+in_uint8 = im2uint8(in);
+color_uint8 = im2uint8(color);
+
+% Initialize the red, green, and blue output channels.
+if ndims(in_uint8) == 2
+    % Input is grayscale.  Initialize all output channels the same.
+    out_red   = in_uint8;
+    out_green = in_uint8;
+    out_blue  = in_uint8;
+else
+    % Input is RGB truecolor.
+    out_red   = in_uint8(:,:,1);
+    out_green = in_uint8(:,:,2);
+    out_blue  = in_uint8(:,:,3);
+end
+
+% Replace output channel values in the mask locations with the appropriate
+% color value.
+out_red(mask)   = color_uint8(1);
+out_green(mask) = color_uint8(2);
+out_blue(mask)  = color_uint8(3);
+
+% Form an RGB truecolor image by concatenating the channel matrices along
+% the third dimension.
+out = cat(3, out_red, out_green, out_blue);

src/main.m

+%% Clear Matlab workspace
+
+clear
+clc
+
+%% User inputs per run
+
+InPath = ['..', filesep, 'dataIn'];
+SavePath = ['..', filesep, 'dataOut'];
+Barcode = {'HO_20180320_Subcloning_snca4correction_1'};
+DisplayChannel = 1;
+
+%% Initialize variables and prepare Java environment
+
+mkdir(SavePath)
+ImPlateSummary = {};
+ObjectsAll = {};
+PreviewCell = {};
+
+%% Analysis
+
+for n = 1:size(Barcode, 2)
+
+    Barcode = Barcode{n};
+    Objects = {};
+
+    %% Document script
+
+    SavePathThisBarcode = [SavePath, filesep, Barcode];
+    AnalysisTimeStamp = datestr(now, 'yyyymmdd_HHMMSS');
+    SavePath = [SavePathThisBarcode, '\Analysis_', AnalysisTimeStamp];
+    mkdir(SavePathThisBarcode)
+    FileNameShort = mfilename;
+    newbackup = sprintf('%s_log.m',[SavePathThisBarcode, filesep, FileNameShort]);
+    FileNameAndLocation = [mfilename('fullpath')];
+    currentfile = strcat(FileNameAndLocation, '.m');
+    copyfile(currentfile,newbackup);
+    Computer = getenv('Computername');
+    Version = ver();
+    save([SavePathThisBarcode, filesep, 'MatlabVersionAndPC.mat'], 'Computer', 'Version')
+    f_LogDependencies(mfilename, SavePathThisBarcode)
+
+    %% Load data
+    [files, info, areaNames, plateLayout, areaMap] = f_Local_Read_Files_Timeseries_Mosaic_BT2('S:\OperaQEHS\OperaArchiveCol', Barcode, 1);
+
+    %% Organize images
+    % find images for a given mosaic
+
+    InfoTable = struct2table(info);
+    IDs = table([1:size(InfoTable, 1)]');
+    IDs.Properties.VariableNames = {'ID'};
+    InfoTable = [InfoTable, IDs];
+    s = 0
+    PlateImChannels = {};
+
+        for row = unique(InfoTable.row)'
+
+            for column = unique(InfoTable.column)'
+
+                if sum((InfoTable.row == row) .* (InfoTable.column == column)) == 0
+                    continue % Replacing try catch
+                end
+
+                s = s + 1
+                ChannelCount = 3; % more channels needed?
+                XYmosaicCells = cell(1,ChannelCount);
+
+                for channel = 1:ChannelCount
+                   [XYmosaic, XYmosaicCell, XPositions, YPositions, ResolutionXY] = f_load_mosaic_image3D_Sparse(files, InfoTable, row, column, channel);
+                   XYmosaicCells{1,channel} = XYmosaicCell;
+                   XYmosaics{1,channel} = XYmosaic;
+                   PlateImChannels{row, column, channel} = imresize(XYmosaics{1, channel}, 0.075);
+                end
+                
+                ch1Green = cell2mat(XYmosaicCells{1,1});% green
+                ch2Red = cell2mat(XYmosaicCells{1,2});% red
+                ch3Blue = cell2mat(XYmosaicCells{1,3});% blue
+                
+                %% Image analysis goes here
+                [ObjectsThisWell, PreviewThisWell] = f_CloneClassifier(row, column, ch1Green, ch2Red, ch3Blue, SavePathThisBarcode, 1/0.075, Barcode);
+                ObjectsAll{row, column} = ObjectsThisWell;
+                PreviewCell{row, column} = PreviewThisWell;
+                
+            end
+        end
+        
+        
+        %% Collect Data
+        
+        Summary = vertcat(ObjectsAll{:});
+        writetable(Summary, [SavePathThisBarcode, filesep, 'summary.csv'])
+        writetable(Summary, [SavePathThisBarcode, filesep, 'summary.xlsx'])
+        
+        %% Plate previews
+        
+        PlateImChannelsEmpyIdx = cellfun(@(x) isempty(x), PlateImChannels);
+        firstBestImSize = size(PlateImChannels{min(find(~PlateImChannelsEmpyIdx))});
+        PlateImChannels(PlateImChannelsEmpyIdx) = {zeros(firstBestImSize, 'uint16')};
+        
+        NotGoodsizedCells = cellfun(@(x) ~min(size(x) == firstBestImSize), PlateImChannels, 'UniformOutput', false); % Remove images of wrong size (file damage exception)
+        PlateImChannels(cell2mat(NotGoodsizedCells)) = {zeros(firstBestImSize, 'uint16')};
+        
+        plateDisplayImCh1  = cell2mat(PlateImChannels(:,:,1));
+        plateDisplayImCh2  = cell2mat(PlateImChannels(:,:,2));
+        plateDisplayImCh3  = cell2mat(PlateImChannels(:,:,3));
+        
+        plateGridMask = zeros(size(plateDisplayImCh1), 'uint16');
+        plateRows = size(PlateImChannels, 1);
+        plateCols = size(PlateImChannels, 2);
+        
+        plateGridMask([1:firstBestImSize(1):(firstBestImSize(1) * plateRows)],:) = 1;
+        plateGridMask(:, [1:firstBestImSize(2):(firstBestImSize(2) * plateCols)],:) = 1;
+        plateGridMask = imdilate(plateGridMask, strel('disk',3));
+        
+        PlatePreviewCh1 = imoverlay(imadjust(plateDisplayImCh1, [0 0.01], [0 1]), plateGridMask, [0 1 1]); % 520/35
+        imwrite(PlatePreviewCh1, [SavePathThisBarcode, filesep, 'ch1_Green.png'])
+        
+        PlatePreviewCh2 = imoverlay(imadjust(plateDisplayImCh2, [0 0.01], [0 1]), plateGridMask, [0 1 1]); % 600/40
+        imwrite(PlatePreviewCh2, [SavePathThisBarcode, filesep, 'ch2_Red.png'])
+        
+        PlatePreviewCh3 = imoverlay(imadjust(plateDisplayImCh3, [0 0.02], [0 1]), plateGridMask, [0 1 1]); % 450/50
+        imwrite(PlatePreviewCh3, [SavePathThisBarcode, filesep, 'ch3_Blue.png'])
+        
+        %% Classification preview
+        
+        firstBestRGBSize = size(PreviewCell{min(find(~PlateImChannelsEmpyIdx))});
+        PreviewCell(max(PlateImChannelsEmpyIdx,[],3)) = {zeros(firstBestRGBSize, 'uint8')};
+        PlatePreviewClassification = cell2mat(PreviewCell);
+        PlatePreviewClassification = imoverlay(PlatePreviewClassification, plateGridMask, [0 1 1]);
+        imwrite(PlatePreviewClassification, [SavePathThisBarcode, filesep, 'Classes.png'])
+        save([SavePathThisBarcode, filesep, 'Workspace.mat'])
+        
+end
+

src/readflexcube.m

+function [ image, meta ] = readflexcube(fname, varargin)
+% READFLEXCUBE reads XYZC data from a Flex file
+% [ cube, meta ] = READFLEXCUBE(fname)
+% * fname is the filename of the flex file
+% * channels is a numeric array with the indexes of the channels to read
+%   Indices start with 1, if not specified read all channels
+% * meta is the BioFormats metadata store object
+
+parser = inputParser;
+addParamValue(parser, 'Info', []);
+addParamValue(parser, 'PlaneCount', []);
+addParamValue(parser, 'ChannelCount', []);
+addParamValue(parser, 'PixelSizeZ', []);
+parse(parser, varargin{:});
+args = parser.Results;
+
+[~, f] = fileattrib(fname);
+fname = f.Name;
+
+if not(isempty(args.Info))
+    info = args.Info;
+else
+    info = imfinfo(fname);
+end
+
+sizeX = info(1).Width;
+sizeY = info(1).Height;
+sizeZ = 1;
+sizeC = 1;
+
+if not(isempty(args.PlaneCount))
+    sizeZ = args.PlaneCount;
+end
+
+if not(isempty(args.ChannelCount))
+    sizeC = args.ChannelCount;
+end
+
+totalcount = numel(info);
+
+if totalcount ~= sizeZ*sizeC
+    if isempty(args.PlaneCount)
+       sizeZ = floor(totalcount / sizeC);
+    elseif isempty(args.ChannelCount)
+       sizeC = floor(totalcount / sizeZ);
+    end
+end
+
+cube = zeros([sizeY(), sizeX(), sizeZ, sizeC], 'uint16');
+
+% getZCTCoords = @(i) [mod(i, sizeZ), ...
+%                      mod(floor(i / sizeZ), sizeC), ...
+%                      floor(floor(i / sizeZ) / sizeC)];
+
+getCZTCoords = @(i) [mod(floor(i / sizeC), sizeZ), ...
+                     mod(i, sizeC), ...
+                     floor(floor(i / sizeC) / sizeZ)];
+
+w = warning('off','all');
+tf = Tiff(fname);
+for i=1:totalcount
+   zct = getCZTCoords(i - 1) + 1;
+   cube(:, :, zct(1), zct(2)) = tf.read();
+   if i < totalcount
+       tf.nextDirectory();
+   end
+end
+warning(w);
+
+image.data = cube;
+image.filename = fname;
+
+xmlstr = info(1).UnknownTags(1).Value;
+xml = xmlreadstring(xmlstr);
+
+image.PixelSizeX = str2double(char(xml.getElementsByTagName('ImageResolutionX').item(0).getTextContent()));
+image.PixelSizeY = str2double(char(xml.getElementsByTagName('ImageResolutionY').item(0).getTextContent()));
+
+if not(isempty(args.PixelSizeZ))
+    pixelZ = args.PixelSizeZ;
+else
+    pixelZ = nan;
+end
+
+image.PixelSizeZ = pixelZ;
+
+meta = xml;
+
+end
+

src/xmlreadstring.m

+function [parseResult,p] = xmlreadstring(stringToParse,varargin)
+%XMLREADSTRING Modified XMLREAD function to read XML data from a string.
+% Author: Luis Cantero.
+% The MathWorks.
+
+p = locGetParser(varargin);
+locSetEntityResolver(p,varargin);
+locSetErrorHandler(p,varargin);
+
+% Parse and return.
+parseStringBuffer = java.io.StringBufferInputStream(stringToParse);
+parseResult = p.parse(parseStringBuffer);
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+function p = locGetParser(args)
+
+p = [];
+for i=1:length(args)
+    if isa(args{i},'javax.xml.parsers.DocumentBuilderFactory')
+        javaMethod('setValidating',args{i},locIsValidating(args));
+        p = javaMethod('newDocumentBuilder',args{i});
+        break;
+    elseif isa(args{i},'javax.xml.parsers.DocumentBuilder')
+        p = args{i};
+        break;
+    end
+end
+
+if isempty(p)
+    parserFactory = javaMethod('newInstance',...
+        'javax.xml.parsers.DocumentBuilderFactory');
+        
+    javaMethod('setValidating',parserFactory,locIsValidating(args));
+    %javaMethod('setIgnoringElementContentWhitespace',parserFactory,1);
+    %ignorable whitespace requires a validating parser and a content model
+    p = javaMethod('newDocumentBuilder',parserFactory);
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+function tf=locIsValidating(args)
+
+tf=any(strcmp(args,'-validating'));
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+function locSetEntityResolver(p,args)
+
+for i=1:length(args)
+    if isa(args{i},'org.xml.sax.EntityResolver')
+        p.setEntityResolver(args{i});
+        break;
+    end
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+function locSetErrorHandler(p,args)
+
+for i=1:length(args)
+    if isa(args{i},'org.xml.sax.ErrorHandler')
+        p.setErrorHandler(args{i});
+        break;
+    end
+end