From 96e0880022ad28a6c2cecdc29e18da86e781194b Mon Sep 17 00:00:00 2001
From: Qianqian Fang <fangqq@gmail.com>
Date: Sat, 21 Dec 2024 22:27:30 -0500
Subject: [PATCH] [feat] add laplacefill.m

---
 fillholes3d.m |   8 ++--
 laplacefill.m | 104 ++++++++++++++++++++++++++++++++++++++++++++++++++
 volgrow.m     |   4 +-
 volshrink.m   |   2 +-
 4 files changed, 111 insertions(+), 7 deletions(-)
 create mode 100644 laplacefill.m

diff --git a/fillholes3d.m b/fillholes3d.m
index 6b1fc7a..1c17818 100644
--- a/fillholes3d.m
+++ b/fillholes3d.m
@@ -21,18 +21,18 @@
 % -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
 %
 
-if (nargin > 1 && numel(maxgap) == 1)
+if (nargin > 1 && numel(maxgap) == 1 && maxgap > 0)
     resimg = volclose(img, maxgap);
 else
     resimg = img;
 end
 
-if (exist('imfill', 'file'))
+if (exist('imfill', 'file') && nargin < 3)
     resimg = imfill(resimg, 'holes');
     return
 end
 
-if (nargin > 1 && numel(maxgap) > 1)
+if (nargin > 1 && ~ischar(maxgap) && numel(maxgap) > 1)
     newimg = zeros(size(resimg) + 2);
 else
     newimg = ones(size(resimg) + 2);
@@ -48,7 +48,7 @@
 end
 
 isseeded = false;
-if (nargin > 1 && numel(maxgap) > 1)
+if (nargin > 1 && ~ischar(maxgap) && numel(maxgap) > 1)
     if (size(maxgap, 2) == 3)
         newimg(sub2ind(size(newimg), maxgap(:, 1) + 1, maxgap(:, 2) + 1, maxgap(:, 3) + 1)) = 1;
     else
diff --git a/laplacefill.m b/laplacefill.m
new file mode 100644
index 0000000..2ff8215
--- /dev/null
+++ b/laplacefill.m
@@ -0,0 +1,104 @@
+function newvol = laplacefill(vol, seeds, solver, tol, maxiter, varargin)
+%
+% newvol = laplacefill(vol, seeds, solver, tol, maxiter, ...)
+%
+% perform a hole-fill or flood-fill in a 2-D or 3-D binary image using
+% Laplace solver
+%
+% author: Qianqian Fang, <q.fang at neu.edu>
+%
+% input:
+%     vol: a 2-D or 3-D binary image
+%     seeds: seeds for flood-fill, 2 or 3-column matrix; if empty, perfrom
+%          hole-fill
+%     solver: (optional) linear solver, if not given, use bicgstab
+%     tol: (optional) linear solver convergence tolerance, default 1e-5
+%     maxiter: (optional) linear solver max iteration, default 3000
+%
+%     additional parameters that are supported by any iterative solver,
+%     such as qmr, pcg, gmres, tfqmr, mldivide, minres etc
+%
+% output:
+%     newvol: the volume image after flood-fill or hole-fill
+%
+% example:
+%     a = zeros(60, 80, 90);
+%     a(20:40, 40:70, 30:80)=1;
+%     a(25:35, 50:60, 50:60)=0;
+%     newvol = laplacefill(a);
+%
+% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
+%
+
+dims = size(vol) + 2;
+
+interiornode = prod(dims - 2);
+
+zeroidx = find(vol == 0);
+
+if (ndims(vol) == 3)
+    [ix, iy, iz] = ndgrid(2:dims(1) - 1, 2:dims(2) - 1, 2:dims(3) - 1);
+else
+    [ix, iy, iz] = ndgrid(2:dims(1) - 1, 2:dims(2) - 1, 2);
+    dims(3) = 3;
+end
+
+ix = ix(zeroidx);
+iy = iy(zeroidx);
+iz = iz(zeroidx);
+
+idx = sub2ind(dims - 2, ix - 1, iy - 1, iz - 1);  % index of all interior nodes in the subset of interior nodes
+idx = idx(:)';
+nnzpos = {idx, idx(ix > 2), idx(ix < dims(1) - 1), idx(iy > 2), idx(iy < dims(2) - 1), idx(iz > 2), idx(iz < dims(3) - 1)};
+
+nonzeroidx = find(vol > 0)';
+
+seedidx = [];
+
+if (nargin > 1 && size(seeds, 2) == 3)
+    seedidx = sub2ind(dims(1:3) - 2, seeds(:, 1), seeds(:, 2), seeds(:, 3));
+end
+
+Amat = sparse([seedidx, nonzeroidx, nnzpos{1}, nnzpos{2} - 1, nnzpos{3} + 1, ...
+               nnzpos{4} - dims(1) + 2, nnzpos{5} + dims(1) - 2, ...
+               nnzpos{6} - (dims(1) - 2) * (dims(2) - 2), nnzpos{7} + (dims(1) - 2) * (dims(2) - 2)], ...
+              [seedidx, nonzeroidx, nnzpos{:}], ...
+              [ones(1, length(seedidx)) ones(1, length(nonzeroidx)) -ones(1, length(idx)) ...
+               (1 / 6) * ones(1, sum(cellfun(@(x) length(x), nnzpos(2:end))))], ...
+              interiornode, interiornode);
+
+if (~isempty(seedidx))
+    b = sparse(seedidx, ones(size(seedidx)), ones(size(seedidx)), interiornode, 1);
+else
+    if (ndims(vol) == 3)
+        boundnode = idx(ix == 2 | iy == 2 | iz == 2 | ix == dims(1) - 1 | iy == dims(2) - 1 | iz == dims(3) - 1);
+    else
+        boundnode = idx(ix == 2 | iy == 2 | ix == dims(1) - 1 | iy == dims(2) - 1);
+    end
+    b = sparse(boundnode, ones(size(boundnode)), -(1 / 6) * ones(size(boundnode)), interiornode, 1);
+end
+
+if (nargin < 4)
+    tol = 1e-10;
+end
+
+if (nargin < 5)
+    maxiter = 3000;
+end
+
+if (nargin > 2)
+    mysolver = str2fun(solver);
+    [newvol, flag] = mysolver(Amat, b, tol, maxiter, varargin{:});
+else
+    [newvol, flag] = bicgstab(Amat, b, tol, maxiter, varargin{:});
+end
+
+if (flag)
+    newvol = gmres(Amat, b, 100, tol, maxiter, varargin{:});
+end
+
+newvol = reshape(full(newvol), size(vol));
+
+if (isempty(seedidx))
+    newvol = ~newvol;
+end
diff --git a/volgrow.m b/volgrow.m
index 43fe736..2350b55 100644
--- a/volgrow.m
+++ b/volgrow.m
@@ -23,7 +23,7 @@
     layer = 1;
 end
 
-if (nargin < 3)
+if (nargin < 3 || isempty(mask))
     if (ndims(vol) == 3)
         mask = zeros(3, 3, 3);
         mask(2, 2, :) = 1;
@@ -37,7 +37,7 @@
 newvol = vol;
 
 for i = 1:layer
-    newvol = (convn(newvol, mask, 'same') > 0);
+    newvol = (convn(single(newvol), single(mask), 'same') > 0);
 end
 
 newvol = double(newvol);
diff --git a/volshrink.m b/volshrink.m
index b2e5ed7..abc1ca8 100644
--- a/volshrink.m
+++ b/volshrink.m
@@ -23,7 +23,7 @@
     layer = 1;
 end
 
-if (nargin < 3)
+if (nargin < 3 || isempty(mask))
     if (ndims(vol) == 3)
         mask = zeros(3, 3, 3);
         mask(2, 2, :) = 1;