Merge pull request #43 from tsipkens/developer

Add draft JOSS paper

Author: tsipkens

+agg/analyze_binary.m

@@ -181,7 +181,7 @@
         %-- Step 3-2: Prepare an image of the isolated aggregate ---------%
         img_binary = zeros(size(img_binary));
         img_binary(CC.PixelIdxList{1,jj}) = 1;
-        Aggs0(jj).binary = logical(img_binary); % store binary image
+        Aggs0(jj).binary = sparse(logical(img_binary)); % store binary image
 
         % Get a cropped version of the aggregate
         % 'autocrop' method included below.
@@ -194,7 +194,7 @@
         img_dilated = imdilate(img_binary,SE);
         img_edge = img_dilated - img_binary;
 
-        [row, col] = find(imcrop(Aggs0(jj).binary, Aggs0(jj).rect));
+        [row, col] = find(imcrop(full(Aggs0(jj).binary), Aggs0(jj).rect));
         Aggs0(jj).length = max((max(row)-min(row)), (max(col)-min(col))) * pixsize(ii);
         Aggs0(jj).width = min((max(row)-min(row)), (max(col)-min(col))) * pixsize(ii);
         Aggs0(jj).aspect_ratio = Aggs0(jj).length / Aggs0(jj).width;
@@ -224,7 +224,11 @@
         [x, y] = find(img_binary ~= 0);
         Aggs0(jj).center_mass = [mean(x); mean(y)];
 
-        if f_plot==1; set(groot,'CurrentFigure',f0); tools.imshow_agg(Aggs0, ii, 0); title(num2str(ii)); drawnow; end
+        if f_plot==1
+            if mod(jj, 10) == 0  % show image after every 10 aggregates processed
+                set(groot,'CurrentFigure',f0); tools.imshow_agg(Aggs0, ii, 0); title(num2str(ii)); drawnow;
+            end
+        end
     end
 
     if f_plot==1; pause(0.05); end % pause very briefly to show overall aggregates
@@ -290,7 +294,6 @@
 y_top = min(max(y)+space,size_img(2)); 
 y_bottom = max(min(y)-space,1);
 
-
 img_binary = img_binary(x_bottom:x_top,y_bottom:y_top);
 img_cropped = img_orig(x_bottom:x_top,y_bottom:y_top);
 rect = [y_bottom,x_bottom,(y_top-y_bottom),(x_top-x_bottom)];

+agg/config/km.v6.1s3.json

@@ -0,0 +1,10 @@
+// VERSION 6.1S
+// AUTHOR: Timothy Sipkens, 2023-04-17
+{
+    "version": "6.1s3",
+    "lvl3": "1:0.002:1.35",  // range of thresholds to build curve
+    "lvl5": 0.1,             // where to threshold on hist curve
+    "lvlfun": "lin",         // type of adjusted threshold
+    "morphsc": 0.5,            // Modify morph. parameter for rolling ball
+    "minsize": 20            // minimum particle size to consider
+}

+agg/config/km.v6.1s4.json

@@ -0,0 +1,10 @@
+// VERSION 6.1S
+// AUTHOR: Timothy Sipkens, 2023-04-17
+{
+    "version": "6.1s4",
+    "lvl3": "1:0.002:1.35",  // range of thresholds to build curve
+    "lvl5": 0.1,             // where to threshold on hist curve
+    "lvlfun": "lin",         // type of adjusted threshold
+    "morphsc": 0.8,            // Modify morph. parameter for rolling ball
+    "minsize": 20            // minimum particle size to consider
+}

+agg/seg_kmeans.m

@@ -94,6 +94,8 @@
 
 
 tools.textheader('k-means');
+disp(['Loaded configuration: ', num2str(opts.version)]);
+disp(' ');
 
 % Loop over images, calling seg function below on each iteration.
 img_binary{n} = []; % pre-allocate cells

+pp/edm_sbs.m

@@ -73,9 +73,9 @@
 
     counts = zeros(length(se_vec),1); % initialize counts
 
-    img_dist = bwdist(~img_binary); % Euclidean distance to outside of aggregate
+    img_dist = bwdist(~full(img_binary)); % Euclidean distance to outside of aggregate
     for ii=1:se_max
-        counts(ii) = nnz(img_dist>se_vec(ii));
+        counts(ii) = nnz(img_dist > se_vec(ii));
             % count the number of non-zero pixels remaining
 
         if counts(ii)==0 % if all of the pixels are gone, exit loop
@@ -102,12 +102,13 @@
     bet = 1.9658; % beta parameter in sigmoid function
     ome = -0.8515; % Omega parameter in sigmoid function
     a = 0.9966;
-    sigmoid = @(x) a./( 1+exp(((log(x(1))-log(dp_bin))./log(x(2))-bet)./ome));
+    sigmoid = @(x) a ./ (1 + ...
+        exp(((log(x(1)) - log(dp_bin)) ./ log(x(2)) - bet) ./ ome));
         % x(1) = dpg, x(2) = spg
 
     opts = optimset('Display','off');
-    x0 = [25,1.5];
-    x1 = lsqnonlin(@(x) (sigmoid(x) - Sa) ./ 100, x0, [], [], opts);
+    x0 = [25, 1.05];
+    x1 = lsqnonlin(@(x) (sigmoid(x) - Sa) ./ 100, x0, [], [100, 1.1], opts);
     x1 = real(x1);
     x1(x1<0) = NaN;
     Sa_fit = sigmoid(x1);

+pp/edm_ws.m

@@ -80,7 +80,7 @@
 
 
     %-- Compute Eucldian distance map ------------------------------------%
-    img_edm = bwdist(~img_binary);
+    img_edm = bwdist(~full(img_binary));
 
     D = max(max(img_edm)) - img_edm;  % EDM to be used for watershed
 

+pp/hough_peng.m

@@ -633,7 +633,7 @@
         cirrad(k) = 0;
     else
         cirrad(k) = SgnCv_LMPos(end);
-        for i_radii = (length(SgnCv_LMPos) - 1) : -1 : 1,
+        for i_radii = (length(SgnCv_LMPos) - 1) : -1 : 1
             circen = [ circen; circen(k,:) ];
             cirrad = [ cirrad; SgnCv_LMPos(i_radii) ];
         end
@@ -643,7 +643,7 @@
 % Output
 varargout{1} = circen;
 varargout{2} = cirrad;
-if nargout > 3,
+if nargout > 3
     varargout{3} = dbg_LMmask;
 end
 

+pp/pcm.m

@@ -65,7 +65,7 @@
 
     %== Step 1: Image preparation ========================================%
     pixsize = Aggs(aa).pixsize; % size of pixels in the image
-    img_binary = imcrop(Aggs(aa).binary, Aggs(aa).rect); % get the binarized image for this aggregate
+    img_binary = imcrop(full(Aggs(aa).binary), Aggs(aa).rect); % get the binarized image for this aggregate
 
     % Get data for this aggregate
     data = Aggs(aa); % initialize data structure for current aggregate

+tools/imshow_agg.m

@@ -81,7 +81,7 @@
             img_binary = or(img_binary,Aggs(aa).binary);
         end
 
-        [~,~,i0] = tools.imshow_binary( ...
+        [~,~,i0] = tools.imshow_binary2( ...
             Aggs(idx1(1)).image, img_binary, opts);
 
         % Make panels bigger.

+tools/imshow_binary.m

@@ -3,69 +3,41 @@
 % Author: Timothy Sipkens, 2019-07-24
 %=========================================================================%
 
-function [h,f,i0] = imshow_binary(imgs, imgs_binary, opts)
+function [h, i1] = imshow_binary(img, img_binary, opts)
 
 %-- Parse inputs ---------------------------------------------------------%
-% Convert images to cells, if they are not already.
-if ~iscell(imgs); imgs = {imgs}; end
-if ~iscell(imgs_binary); imgs_binary = {imgs_binary}; end
-
-if length(imgs)>24 % only plot up to 24 images
-    imgs = imgs(1:24);
-    imgs_binary = imgs_binary(1:24);
-end
-n_imgs = length(imgs); % number of images
+if iscell(img); img = img{1}; end  % select only first image
+if iscell(img_binary); img_binary = img_binary{1}; end
 
 if ~exist('opts','var'); opts = struct(); end
 if ~isfield(opts,'cmap')
-    opts.cmap = ones(max(max(imgs_binary{1})),1) * [0.12,0.59,0.96];
+    opts.cmap = ones(max(max(img_binary)),1) * [0,1,1];
         % blue, repeated for as many labels as exist in img_binary
 end
 if ~isfield(opts,'f_outline'); opts.f_outline = 1; end
 if ~isfield(opts,'label_alpha'); opts.label_alpha = 0.7; end
 %-------------------------------------------------------------------------%
 
 
+t0 = labeloverlay(img, img_binary, ...
+    'Transparency', opts.label_alpha,...
+    'Colormap', opts.cmap);
 
-% If more than one image, prepare to tile and maximize figure.
-if n_imgs>1
-    N1 = floor(sqrt(n_imgs));
-    N2 = ceil(n_imgs/N1);
-    subplot(N1, N2, 1);
-end
-
-
-for ii=1:n_imgs % loop over images
-    
-    if n_imgs>1
-        subplot(N1, N2, ii);
-        title(num2str(ii));
-    end
-    
-    t0 = labeloverlay(imgs{ii}, imgs_binary{ii}, ...
-        'Transparency', opts.label_alpha,...
-        'Colormap', opts.cmap);
-    
-    if ~opts.f_outline
-        i1 = t0;
-    else % if adding an outline
-        img_edge = edge(imgs_binary{ii},'sobel');
-        se = strel('disk',1);
-        img_dilated = imdilate(img_edge,se);
-            % use dilation to strengthen the aggregate's outline
-        i1 = uint8(~img_dilated).*t0;
-            % adds borders to labeled regions
-    end
-    
-    tools.imshow(i1); % show labelled image
-    
-    if ii==1; i0 = i1; end % store first figure panel for output
+if ~opts.f_outline
+    i1 = t0;
 
+else % if adding an outline
+    img_edge = edge(img_binary,'sobel');
+    se = strel('disk',1);
+    img_dilated = imdilate(img_edge,se);
+        % use dilation to strengthen the aggregate's outline
+    i1 = uint8(~img_dilated).*t0;
+        % adds borders to labeled regions
 end
 
+tools.imshow(i1); % show labelled image
 
-if nargout>0; h = gca; end % organize outputs
-if nargout>1; f = gcf; end
+if nargout>0; h = gca; end % organize output
 
 end
 

+tools/imshow_binary2.m

@@ -0,0 +1,49 @@
+
+% INSHOW_BINARY  Plot original image with binary mask overlayed
+% Author: Timothy Sipkens, 2019-07-24
+%=========================================================================%
+
+function [h,f,i0] = imshow_binary2(imgs, imgs_binary, varargin)
+
+%-- Parse inputs ---------------------------------------------------------%
+% Convert images to cells, if they are not already.
+if ~iscell(imgs); imgs = {imgs}; end
+if ~iscell(imgs_binary); imgs_binary = {imgs_binary}; end
+
+if length(imgs)>24 % only plot up to 24 images
+    imgs = imgs(1:24);
+    imgs_binary = imgs_binary(1:24);
+end
+n_imgs = length(imgs); % number of images
+%-------------------------------------------------------------------------%
+
+
+
+% If more than one image, prepare to tile and maximize figure.
+if n_imgs>1
+    N1 = floor(sqrt(n_imgs));
+    N2 = ceil(n_imgs/N1);
+    subplot(N1, N2, 1);
+end
+
+
+for ii=1:n_imgs % loop over images
+    
+    if n_imgs>1
+        subplot(N1, N2, ii);
+        title(num2str(ii));
+    end
+    
+    [~, i1] = tools.imshow_binary(imgs{ii}, imgs_binary{ii}, varargin{:});
+    
+    if ii==1
+        if nargout>0; h = gca; end  % organize outputs
+        i0 = i1;  % store first figure panel for output
+    end
+    
+end
+
+if nargout>1; f = gcf; end
+
+end
+

+tools/load_imgs.m

@@ -291,10 +291,14 @@
         bw1 = im2bw(1 - double(Imgs(jj).raw) ./ ...
             max(max(double(Imgs(jj).raw))), 0.98);
 
-        % If no text, loop for white text isntead of black.
-        if nnz(bw1) < 500
+        % If no text, loop for white text instead of black.
+        if nnz(bw1) / numel(bw1) < 0.01
             bw1 = ~im2bw(1 - double(Imgs(jj).raw) ./ ...
-                max(max(double(Imgs(jj).raw))), 0.01);
+                max(max(double(Imgs(jj).raw))), 0.08);
+
+            se = strel('disk', 1);
+            bw1 = imclose(bw1, se);
+            bw1 = imopen(bw1, se);
         end
 
         % Remove any small regions below a certain number of pixels.

.github/workflows/draft-pdf.yml

@@ -0,0 +1,19 @@
+on: [push]
+
+jobs:
+  paper:
+    runs-on: ubuntu-latest
+    name: Paper Draft
+    steps:
+      - name: Checkout
+        uses: actions/checkout@v3
+      - name: Build draft PDF
+        uses: openjournals/openjournals-draft-action@master
+        with:
+          journal: joss
+          paper-path: docs/paper/paper.md
+      - name: Upload
+        uses: actions/upload-artifact@v1
+        with:
+          name: paper
+          path: docs/paper/paper.pdf

.gitignore

@@ -9,8 +9,8 @@ results/
 results\
 +working/
 +\working
-figs/
-figs\
+docs/figs/
+docs\figs\
 
 temp/
 temp\
@@ -173,3 +173,4 @@ cython_debug/
 
 carboseg/FPN-resnet50-imagenet.onnx
 +agg/config/km.v6.1s2.json
+docs/paper/instructions.txt

README.md

@@ -393,7 +393,7 @@ This program contains very significantly modified versions of the code distribut
 
 Also included with this program is the MATLAB code of [Kook et al. (2015)][kook], modified to accommodate the expected inputs and outputs common to the other functions.
 
-This code also contain an adaptation of **EDM-SBS** method of [Bescond et al. (2014)][bescond]. We thank the authors, in particular Jérôme Yon, for their help in understanding their original [Scilab code and ImageJ plugin](http://www.coria.fr/spip.php?article910). Modifications to allow the method to work directly on binary images (rather than a custom output from ImageJ) and to integrate the method into the MATLAB environment may present some minor compatibility issues, but allows use of the aggregate segmentation methods given in the **agg** package. 
+This code also contain an adaptation of **EDM-SBS** method of [Bescond et al. (2014)][bescond]. We thank the authors, in particular Jérôme Yon, for their help in understanding their original [Scilab code and ImageJ plugin](https://www.coria.fr/en/edm-sbs-automated-analysis-of-tem-images/). Modifications to allow the method to work directly on binary images (rather than a custom output from ImageJ) and to integrate the method into the MATLAB environment may present some minor compatibility issues, but allows use of the aggregate segmentation methods given in the **agg** package. 
 
 The **carboseg** method follows from collaborative work with [Max Frei](https://github.com/maxfrei750) and is associated with [Sipkens et al. (2021)][ptech.cnn].