Mercurial > hg > octave-image
changeset 878:6ef0d8db565a
maint: Periodic merge of stable to default.
| author | Carnë Draug <carandraug@octave.org> |
|---|---|
| date | Sat, 08 Mar 2014 01:57:34 +0000 |
| parents | f1cf02800a87 (current diff) 72930d024906 (diff) |
| children | 090388cdb583 |
| files | NEWS |
| diffstat | 8 files changed, 467 insertions(+), 175 deletions(-) [+] |
line wrap: on
line diff
--- a/.hgtags +++ b/.hgtags @@ -32,3 +32,4 @@ 587e29845b1d0b82d80ac460a832b40053150e76 Octave-Forge-2001.11.02 dbe605c10a08046e56b025926d1c3ec08136d952 snapshot-2.1.1 b0977102e7f65a2b220b07e46b64486b5ccfbc35 release-2.2.0 +fabc22044c6b7f1c19af0f1942cfc4c18e90280e release-2.2.1
--- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Name: Image -Version: 2.2.0 -Date: 2014-01-08 +Version: 2.2.1 +Date: 2014-03-08 Author: various authors Maintainer: Carnë Draug <carandraug+dev@gmail.com> Title: Image Processing
--- a/NEWS +++ b/NEWS @@ -6,9 +6,6 @@ weighted by how much of them is covered by the disk. Note that this change is backwards incompatible. - ** The 'disk' filter of the function fspecial has been changed in a - backwards incompatible way. - ** The following functions will display the output image as a figure instead of printing to the command line, when there are no output arguments: @@ -24,6 +21,25 @@ impad iptcheckstrs imrotate_Fourier uintlut + + Summary of important user-visible changes for image 2.2.1 (2014/03/08): +------------------------------------------------------------------------- + + ** imcrop had many alternative interfaces added for more flexibility. + Added support in the input for indexed images, figures handles, + N-dimensional images, and specific bounding box vector for a + non-interactive usage. Output can now also return the bounding box used + for the cropping in addition to the cropped image. It will no longer + loop forever until it gets two valid coordinates for the bounding box. + + ** Fixed bug in imcomplement to compute the complement of signed integers + correctly. + + ** Fix imrotate to handle RGB images. + + ** Fix regression in bwdist when calculating the closest pixel map. + + Summary of important user-visible changes for image 2.2.0 (2014/01/08): -------------------------------------------------------------------------
--- a/inst/imcomplement.m +++ b/inst/imcomplement.m @@ -1,4 +1,5 @@ ## Copyright (C) 2008 Søren Hauberg <soren@hauberg.org> +## Copyright (C) 2014 Carnë Draug <carandraug@octave.org> ## ## This program is free software; you can redistribute it and/or modify it under ## the terms of the GNU General Public License as published by the Free Software @@ -14,33 +15,73 @@ ## this program; if not, see <http://www.gnu.org/licenses/>. ## -*- texinfo -*- -## @deftypefn {Function File} @var{B} = imcomplement(@var{A}) -## Computes the complement image. Intuitively this corresponds to the intensity -## of bright and dark regions being reversed. +## @deftypefn {Function File} {} imcomplement (@var{A}) +## Compute image complement or negative. +## +## Intuitively this corresponds to the intensity of bright and dark +## regions being reversed. The exact operation performed is dependent +## on the class of the image. ## -## For binary images, the complement is computed as @code{!@var{A}}, for floating -## point images it is computed as @code{1 - @var{A}}, and for integer images as -## @code{intmax(class(@var{A})) - @var{A}}. +## @table @asis +## @item floating point +## Since floating point images are meant to have values in the range [0 1], +## this is equivalent @code{A -1}. This leads to values within the range +## to be inverted while others to equally distance from the limits. +## +## @item logical +## Equivalent to @code{! A} +## +## @item integer +## Inverts the values within the range of the data. This is a different +## depending whether it's signed or unsigned integers class. +## @end table +## ## @seealso{imadd, imdivide, imlincomb, immultiply, imsubtract} ## @end deftypefn -function B = imcomplement(A) - ## Check input +function B = imcomplement (A) + if (nargin != 1) - error("imcomplement: not enough input arguments"); - endif - if (!ismatrix(A)) - error("imcomplement: input must be an array"); + print_usage (); endif - ## Take action depending on the class of A - if (isa(A, "double") || isa(A, "single")) + if (isfloat (A)) B = 1 - A; - elseif (islogical(A)) - B = !A; - elseif (isinteger(A)) - B = intmax(class(A)) - A; + elseif (islogical (A)) + B = ! A; + elseif (any (isa (A, {"int8", "int16", "int32", "int64"}))) + ## Signed integers are different. We have all this extra work to avoid + ## conversion into another class which may not fit in memory. And seems + ## like it it still performs faster. + cl = class (A); + B = zeros (size (A), cl); + m = (A != intmin (cl)); + + B(! m) = intmax (cl); + B(m) = -A -1; + + elseif (isinteger (A)) + B = intmax (class (A)) - A; else - error("imcomplement: unsupported input class: '%s'", class(A)); + error("imcomplement: A must be an image but is of class `%s'", class (A)); endif + endfunction + +%!assert (imcomplement (10), -9); +%!assert (imcomplement (single (10)), single (-9)); +%!assert (imcomplement (0.2), 0.8); +%!assert (imcomplement (uint8 (0)), uint8 (255)); +%!assert (imcomplement (uint8 (1)), uint8 (254)); +%!assert (imcomplement (uint16 (0)), uint16 (65535)); +%!assert (imcomplement (uint16 (1)), uint16 (65534)); + +%!assert (imcomplement (int8 (-128)), int8 ( 127)); +%!assert (imcomplement (int8 ( 127)), int8 (-128)); +%!assert (imcomplement (int16 (-1)), int16 ( 0)); +%!assert (imcomplement (int16 ( 0)), int16 (-1)); +%!assert (imcomplement (int16 ( 1)), int16 (-2)); + +%!assert (imcomplement ([true false true]), [false true false]) + +%!error <must be an image> imcomplement ("not an image")
--- a/inst/imcrop.m +++ b/inst/imcrop.m @@ -1,5 +1,4 @@ -## Copyright (C) 2012 Pablo Rossi <prossi@ing.unrc.edu.ar> -## Copyright (C) 2012 Carnë Draug <carandraug+dev@gmail.com> +## Copyright (C) 2014 Carnë Draug <carandraug+dev@gmail.com> ## ## This program is free software; you can redistribute it and/or modify it under ## the terms of the GNU General Public License as published by the Free Software @@ -15,36 +14,211 @@ ## this program; if not, see <http://www.gnu.org/licenses/>. ## -*- texinfo -*- -## @deftypefn {Function File} @var{cropped} = imcrop (@var{Img}) +## @deftypefn {Function File} {} imcrop () +## @deftypefnx {Function File} {} imcrop (@var{img}) +## @deftypefnx {Function File} {} imcrop (@var{ind}, @var{cmap}) +## @deftypefnx {Function File} {} imcrop (@var{h}) +## @deftypefnx {Function File} {} imcrop (@dots{}, @var{rect}) +## @deftypefnx {Function File} {[@var{cropped}] =} imcrop (@dots{}) +## @deftypefnx {Function File} {[@var{cropped}, @var{rect}] =} imcrop (@dots{}) +## @deftypefnx {Function File} {[@var{x}, @var{y}, @var{cropped}, @var{rect}] =} imcrop (@dots{}) ## Crop image. ## -## Displays the image @var{Img} in a figure window and waits for user to select -## two points to define the bounding box. First click on the top left and then -## on the bottom right corner of the region. The function will not return until -## two valid points in the correct order are selected. +## Displays the image @var{img} in a figure window and waits for the user to +## select two points defining a bounding box. First click on the top left +## and then on the bottom right corner of the region. For an indexed image, a +## corresponding colormap can be specified in @var{cmap}. For multi-dimensional +## images (each 2D image is concatenated in the 4th dimension), only the +## first image is displayed. +## +## if no image data is given, uses the current figure or figure from graphics +## handle @var{h}. +## +## Non-interactive usage is supported if the last input argument is 4 element +## vector @var{rect} defining the region of interest. The first two elements +## specify the initial @var{x_ini} and @var{y_ini} coordinates, and the last +## two the @var{width} and @var{height}, i.e., +## @code{@var{rect} = [@var{x_ini} @var{y_ini} @var{width} @var{height}]}. +## Note how this the opposite of the majority of Octave indexing rules where +## rows come before columns. +## +## Returns the @var{cropped} image and a vector @var{rect} with the +## coordinates and size for @var{cropped}. If more than 3 output arguments +## are requested, also returns the @var{x} and @var{y} data that define +## the coordinate system. ## -## Returns the @var{cropped} image. +## @emph{Note}: the values in @var{rect} are not necessarily integer values +## and can't always be used directly as index values for other images. To +## crop the same region from a multiple images of the same size, either using +## a multi-dimensional image: ## -## @seealso{imshow} +## @example +## @group +## nd_img = cat (4, img1, img2, img3, img4); +## croped = imcrop (nd_img); +## cropped_1 = cropped(:,:,:,1); +## cropped_2 = cropped(:,:,:,2); +## cropped_3 = cropped(:,:,:,3); +## cropped_4 = cropped(:,:,:,4); +## @end group +## @end example +## +## or multiple calls to @code{imcrop}: +## +## @example +## @group +## [croped_1, rect] = imcrop (img1); +## cropped_2 = imcrop (img2, rect); +## cropped_3 = imcrop (img3, rect); +## cropped_4 = imcrop (img4, rect); +## @end group +## @end example +## +## @seealso{impixel, imshow} ## @end deftypefn -function col = imcrop (Img) +## TODO not yet implemented +## @deftypefnx {Function File} {} imcrop (@var{xData}, @var{yData}, @dots{}) + +function varargout = imcrop (varargin) + + ## Screw Matlab and their over complicated API's! How can we properly + ## parse all the possible alternative calls? See + ## http://www.youtube.com/watch?v=1oZWacjmYm8 to understand how such + ## API's develop. + + ## There is no check for this things, anything is valid. We (Octave) + ## are at least checking the number of elements otherwise the input + ## parsing would be horrible. + valid_rect = @(x) numel (x) == 4; + valid_system = @(x) numel (x) == 2; + + rect = []; + interactive = true; # is interactive usage + alt_system = false; # was an alternative coordinate system requested? + from_fig = false; # do we have the image data or need to fetch from figure? + + if (nargin > 5) + print_usage (); + endif + + rect = []; + if (numel (varargin) > 1 && valid_rect (varargin{end})) + interactive = false; + rect = varargin{end}; + varargin(end) = []; + endif - handle = imshow (Img); - [a, b] = size (Img); + xdata = []; + ydata = []; + if (numel (varargin) > 2 && valid_system (varargin{1}) && valid_system (varargin{2})) + ## requested messy stuff + ## we should probably reuse part of what impixel does + alt_system = true; + xdata = varargin{1}; + ydata = varargin{2}; + varargin([1 2]) = []; + error ("imcrop: messing around with coordinate system is not implemented"); + endif + + ## After we remove all that extra stuff, we are left with the image + fnargin = numel (varargin); + if (fnargin > 2) + print_usage (); + elseif (fnargin == 0) + ## use current figure + from_fig = true; + h = gcf (); + elseif (fnargin == 1 && ishandle (varargin{1})) + ## use specified figure + from_fig = true; + h = varargin{1}; + elseif (interactive) + ## leave input check to imshow + h = nd_imshow (varargin{:}); + elseif (isimage (varargin{1})) + ## we have the image data and it's not interactive, so there is + ## nothing to do. We only check the minimum in the image. + else + print_usage (); + endif - do - [hl, rd] = ginput(2); - if (hl(1) <= 1), hl(1) = 1; endif - if (rd(1) <= 1), rd(1) = 1; endif - if (hl(2) >= b), hl(2) = b; endif - if (rd(2) >= a), rd(2) = a; endif - until (hl(1) < hl(2) || rd(1) < rd(2)) - ## should we close the image after? Anyway, close does not accept the handle - ## since the handle from imshow is not a figure handle -# close (handle); + if (from_fig) + cdata = get (h, "cdata"); + xdata = get (h, "xdata"); + ydata = get (h, "ydata"); + else + cdata = varargin{1}; + if (! alt_system) + xdata = [1 columns(cdata)]; + ydata = [1 rows(cdata)]; + endif + endif + + ## Finally, crop the image + if (interactive) + [x, y] = ginput (2); + rect = [x(1) y(1) x(2)-x(1) y(2)-y(1)]; + endif + i_ini = round ([rect(1) rect(2)]); + i_end = round ([rect(1)+rect(3) rect(2)+rect(4)]); + img = cdata(i_ini(2):i_end(2), i_ini(1):i_end(1),:,:); # don't forget RGB and ND images + + ## Even the API for the output is complicated + if (nargout == 0 && interactive) + figure (); + ## In case we have a colormap or something like that, use + ## it again when displaying the cropped image. + nd_imshow (img, varargin{2:end}); + elseif (nargout < 3) + varargout{1} = img; + varargout{2} = rect; + else + varargout{1} = xdata; + varargout{2} = ydata; + varargout{3} = img; + varargout{4} = rect; + endif + +endfunction - hl = floor (hl); - rd = floor (rd); - col = Img(rd(1):rd(2), hl(1):hl(2)); +## shadows core function to support ND image. If we have one, use +## the first frame only +function h = nd_imshow (varargin) + size (varargin{1}); + h = imshow (varargin{1}(:,:,:,1), varargin{2:end}); endfunction + +## test typical non-interactive usage, grayscale image +%!test +%! a = randi (255, [100 100]); +%! rect = [20 30 3 5]; +%! assert (nthargout ([1 2], @imcrop, a, rect), {a(30:35, 20:23) rect}); +%! assert (nthargout (2, @imcrop, a, rect), rect); +%! assert (nthargout ([3 4], 4, @imcrop, a, rect), {a(30:35, 20:23) rect}); + +## test typical non-interactive usage, RGB image +%!test +%! rgb = randi (255, [100 100 3]); +%! rect = [20 30 3 5]; +%! assert (nthargout ([1 2], @imcrop, rgb, rect), {rgb(30:35, 20:23,:) rect}); +%! assert (nthargout (2, @imcrop, rgb, rect), rect); +%! assert (nthargout ([3 4], 4, @imcrop, rgb, rect), {rgb(30:35, 20:23,:) rect}); + +## test typical non-interactive usage, indexed image +%!test +%! a = randi (255, [100 100]); +%! rect = [20 30 3 5]; +%! cmap = jet (255); +%! assert (nthargout ([1 2], @imcrop, a, cmap, rect), {a(30:35, 20:23) rect}); +%! assert (nthargout (2, @imcrop, a, cmap, rect), rect); +%! assert (nthargout ([3 4], 4, @imcrop, a, cmap, rect), {a(30:35, 20:23) rect}); + +## test typical non-interactive usage, logical image +%!test +%! a = rand (100) > 0.5; +%! rect = [20 30 3 5]; +%! assert (nthargout ([1 2], @imcrop, a, rect), {a(30:35, 20:23) rect}); +%! assert (nthargout (2, @imcrop, a, rect), rect); +%! assert (nthargout ([3 4], 4, @imcrop, a, rect), {a(30:35, 20:23) rect}); +
--- a/inst/imrotate.m +++ b/inst/imrotate.m @@ -61,7 +61,7 @@ if (nargin < 2 || nargin > 5) print_usage (); - elseif (! isimage (imgPre) || ndims (imgPre) != 2) + elseif (! isimage (imgPre)) error ("imrotate: IMGPRE must be a grayscale or RGB image.") elseif (! isscalar (thetaDeg)) error("imrotate: THETA must be a scalar"); @@ -118,12 +118,12 @@ nRot90 = mod(thetaDeg, 360) / 90; if (mod(thetaDeg, 180) == 0 || sizePre(1) == sizePre(2) || strcmpi(bbox, "loose")) - imgPost = rot90(imgPre, nRot90); + imgPost = rotdim (imgPre, nRot90, [1 2]); return; elseif (mod(sizePre(1), 2) == mod(sizePre(2), 2)) ## Here, bbox is "crop" and the rotation angle is +/- 90 degrees. ## This works only if the image dimensions are of equal parity. - imgRot = rot90(imgPre, nRot90); + imgRot = rotdim (imgPre, nRot90, [1 2]); imgPost = zeros(sizePre); hw = min(sizePre) / 2 - 0.5; imgPost (round(oPost(2) - hw) : round(oPost(2) + hw), @@ -139,7 +139,7 @@ ## caller who wants to avoid this should ensure that the image ## dimensions are of equal parity. endif - end + endif ## Now the actual rotations happen if (strcmpi (interp, "fourier")) @@ -188,14 +188,14 @@ phi = theta; elseif ( theta>45 && theta<=135 ) phi = theta - 90; - f = rot90(f,1); + f = rotdim(f,1, [1 2]); perp = 1; elseif ( theta>135 && theta<=225 ) phi = theta - 180; - f = rot90(f,2); + f = rotdim(f,2, [1 2]); elseif ( theta>225 && theta<=315 ) phi = theta - 270; - f = rot90(f,3); + f = rotdim(f,3, [1 2]); perp = 1; else phi = theta; @@ -293,57 +293,82 @@ endfunction -%!test -%! ## Verify minimal loss across six rotations that add up to 360 +/- 1 deg.: -%! methods = { "nearest", "bilinear", "bicubic", "Fourier" }; -%! angles = [ 59 60 61 ]; -%! tolerances = [ 7.4 8.5 8.6 # nearest -%! 3.5 3.1 3.5 # bilinear -%! 2.7 2.0 2.7 # bicubic -%! 2.7 1.6 2.8 ]/8; # Fourier -%! -%! # This is peaks(50) without the dependency on the plot package -%! x = y = linspace(-3,3,50); -%! [X,Y] = meshgrid(x,y); -%! x = 3*(1-X).^2.*exp(-X.^2 - (Y+1).^2) ... -%! - 10*(X/5 - X.^3 - Y.^5).*exp(-X.^2-Y.^2) ... -%! - 1/3*exp(-(X+1).^2 - Y.^2); -%! -%! x -= min(x(:)); # Fourier does not handle neg. values well -%! x = x./max(x(:)); -%! for m = 1:(length(methods)) -%! y = x; -%! for i = 1:5 -%! y = imrotate(y, 60, methods{m}, "crop", 0); -%! end -%! for a = 1:(length(angles)) -%! assert(norm((x - imrotate(y, angles(a), methods{m}, "crop", 0)) -%! (10:40, 10:40)) < tolerances(m,a)); -%! end -%! end +#%!test +#%! ## Verify minimal loss across six rotations that add up to 360 +/- 1 deg.: +#%! methods = { "nearest", "bilinear", "bicubic", "Fourier" }; +#%! angles = [ 59 60 61 ]; +#%! tolerances = [ 7.4 8.5 8.6 # nearest +#%! 3.5 3.1 3.5 # bilinear +#%! 2.7 2.0 2.7 # bicubic +#%! 2.7 1.6 2.8 ]/8; # Fourier +#%! +#%! # This is peaks(50) without the dependency on the plot package +#%! x = y = linspace(-3,3,50); +#%! [X,Y] = meshgrid(x,y); +#%! x = 3*(1-X).^2.*exp(-X.^2 - (Y+1).^2) ... +#%! - 10*(X/5 - X.^3 - Y.^5).*exp(-X.^2-Y.^2) ... +#%! - 1/3*exp(-(X+1).^2 - Y.^2); +#%! +#%! x -= min(x(:)); # Fourier does not handle neg. values well +#%! x = x./max(x(:)); +#%! for m = 1:(length(methods)) +#%! y = x; +#%! for i = 1:5 +#%! y = imrotate(y, 60, methods{m}, "crop", 0); +#%! end +#%! for a = 1:(length(angles)) +#%! assert(norm((x - imrotate(y, angles(a), methods{m}, "crop", 0)) +#%! (10:40, 10:40)) < tolerances(m,a)); +#%! end +#%! end -%!xtest -%! ## Verify exactness of near-90 and 90-degree rotations: -%! X = rand(99); -%! for angle = [90 180 270] -%! for da = [-0.1 0.1] -%! Y = imrotate(X, angle + da , "nearest", :, 0); -%! Z = imrotate(Y, -(angle + da), "nearest", :, 0); -%! assert(norm(X - Z) == 0); # exact zero-sum rotation -%! assert(norm(Y - imrotate(X, angle, "nearest", :, 0)) == 0); # near zero-sum -%! end -%! end +#%!xtest +#%! ## Verify exactness of near-90 and 90-degree rotations: +#%! X = rand(99); +#%! for angle = [90 180 270] +#%! for da = [-0.1 0.1] +#%! Y = imrotate(X, angle + da , "nearest", :, 0); +#%! Z = imrotate(Y, -(angle + da), "nearest", :, 0); +#%! assert(norm(X - Z) == 0); # exact zero-sum rotation +#%! assert(norm(Y - imrotate(X, angle, "nearest", :, 0)) == 0); # near zero-sum +#%! end +#%! end + +#%!test +#%! ## Verify preserved pixel density: +#%! methods = { "nearest", "bilinear", "bicubic", "Fourier" }; +#%! ## This test does not seem to do justice to the Fourier method...: +#%! tolerances = [ 4 2.2 2.0 209 ]; +#%! range = 3:9:100; +#%! for m = 1:(length(methods)) +#%! t = []; +#%! for n = range +#%! t(end + 1) = sum(imrotate(eye(n), 20, methods{m}, :, 0)(:)); +#%! end +#%! assert(t, range, tolerances(m)); +#%! end %!test -%! ## Verify preserved pixel density: -%! methods = { "nearest", "bilinear", "bicubic", "Fourier" }; -%! ## This test does not seem to do justice to the Fourier method...: -%! tolerances = [ 4 2.2 2.0 209 ]; -%! range = 3:9:100; -%! for m = 1:(length(methods)) -%! t = []; -%! for n = range -%! t(end + 1) = sum(imrotate(eye(n), 20, methods{m}, :, 0)(:)); -%! end -%! assert(t, range, tolerances(m)); -%! end +%! a = reshape (1:18, [2 3 3]); +%! +%! a90(:,:,1) = [5 6; 3 4; 1 2]; +%! a90(:,:,2) = a90(:,:,1) + 6; +%! a90(:,:,3) = a90(:,:,2) + 6; +%! +%! a180(:,:,1) = [6 4 2; 5 3 1]; +%! a180(:,:,2) = a180(:,:,1) + 6; +%! a180(:,:,3) = a180(:,:,2) + 6; +%! +%! am90(:,:,1) = [2 1; 4 3; 6 5]; +%! am90(:,:,2) = am90(:,:,1) + 6; +%! am90(:,:,3) = am90(:,:,2) + 6; +%! +%! assert (imrotate (a, 0), a); +%! assert (imrotate (a, 90), a90); +%! assert (imrotate (a, -90), am90); +%! assert (imrotate (a, 180), a180); +%! assert (imrotate (a, -180), a180); +%! assert (imrotate (a, 270), am90); +%! assert (imrotate (a, -270), a90); +%! assert (imrotate (a, 360), a); +
--- a/src/bwdist.cc +++ b/src/bwdist.cc @@ -41,7 +41,7 @@ */ void edtfunc (float (*func)(short int, short int), - const Matrix &img, + const boolNDArray &img, short *distx, short *disty) { @@ -49,45 +49,47 @@ const int h = img.rows (); const int numel = img.numel (); - int x, y, i; + // Initialize the distance images to be all large values + const bool* elem = img.fortran_vec (); + for (octave_idx_type i = 0; i < numel; i++) + { + if(elem[i]) + { + distx[i] = 0; + disty[i] = 0; + } + else + { + distx[i] = 32000; // Large but still representable in a short, and + disty[i] = 32000; // 32000^2 + 32000^2 does not overflow an int + } + } + float olddist2, newdist2, newdistx, newdisty; bool changed; // Initialize index offsets for the current image width - const int offset_u = -w; - const int offset_ur = -w+1; + const int offset_u = -h; + const int offset_ur = -h+1; const int offset_r = 1; - const int offset_rd = w+1; - const int offset_d = w; - const int offset_dl = w-1; + const int offset_rd = h+1; + const int offset_d = h; + const int offset_dl = h-1; const int offset_l = -1; - const int offset_lu = -w-1; - - // Initialize the distance images to be all large values - for (i = 0; i < numel; i++) - { - if(img(i) == 0.0) - { - distx[i] = 32000; // Large but still representable in a short, and - disty[i] = 32000; // 32000^2 + 32000^2 does not overflow an int - } - else - { - distx[i] = 0; - disty[i] = 0; - } - } + const int offset_lu = -h-1; // Perform the transformation + int x, y, i; do { changed = false; // Scan rows, except first row - for (y = 1; y < h; y++) + for (y = 1; y < w; y++) { + OCTAVE_QUIT; // move index to leftmost pixel of current row - i = y*w; + octave_idx_type i = y*h; /* scan right, propagate distances from above & left */ @@ -119,9 +121,8 @@ i++; /* Middle pixels have all neighbors */ - for(x=1; x<w-1; x++, i++) + for(x=1; x<h-1; x++, i++) { - OCTAVE_QUIT; olddist2 = (*func)(distx[i], disty[i]); if(olddist2 == 0) continue; // Already zero distance @@ -209,12 +210,11 @@ /* Move index to second rightmost pixel of current row. */ /* Rightmost pixel is skipped, it has no right neighbor. */ - i = y*w + w-2; + i = y*h + h-2; /* scan left, propagate distance from right */ - for(x=w-2; x>=0; x--, i--) + for(x=h-2; x>=0; x--, i--) { - OCTAVE_QUIT; olddist2 = (*func)(distx[i], disty[i]); if(olddist2 == 0) continue; // Already zero distance @@ -231,11 +231,11 @@ } /* Scan rows in reverse order, except last row */ - for(y=h-2; y>=0; y--) + for(y=w-2; y>=0; y--) { OCTAVE_QUIT; /* move index to rightmost pixel of current row */ - i = y*w + w-1; + i = y*h + h-1; /* Scan left, propagate distances from below & right */ @@ -267,9 +267,8 @@ i--; /* Middle pixels have all neighbors */ - for(x=w-2; x>0; x--, i--) + for(x=h-2; x>0; x--, i--) { - OCTAVE_QUIT; olddist2 = (*func)(distx[i], disty[i]); if(olddist2 == 0) continue; // Already zero distance @@ -356,10 +355,9 @@ /* Move index to second leftmost pixel of current row. */ /* Leftmost pixel is skipped, it has no left neighbor. */ - i = y*w + 1; - for(x=1; x<w; x++, i++) + i = y*h + 1; + for(x=1; x<h; x++, i++) { - OCTAVE_QUIT; /* scan right, propagate distance from left */ olddist2 = (*func)(distx[i], disty[i]); if(olddist2 == 0) continue; // Already zero distance @@ -382,7 +380,8 @@ // The different functions used to calculate distance, as a // class so its typename can be used for edtfunc template -static float euclidean (short x, short y) +static float +euclidean (short x, short y) { // the actual euclidean distance, is the square root of this. But // squaring does not change the order of the distances, so we can @@ -390,40 +389,48 @@ return ((int)(x)*(x) + (y)*(y)); } -static float chessboard (short x, short y) +static float +chessboard (short x, short y) { return std::max (abs (y), abs (x)); } -static float cityblock (short x, short y) +static float +cityblock (short x, short y) { return abs (x) + abs (y); } -static float quasi_euclidean (short x, short y) +static float +quasi_euclidean (short x, short y) { static const float sqrt2_1 = sqrt (2) - 1; return abs(x)>abs(y) ? (abs(x) + sqrt2_1 * abs(y)) : (sqrt2_1 * abs(x) + abs(y)) ; } -FloatMatrix calc_distances (float (*func)(short, short), - Matrix bw, - short *xdist, - short *ydist) +static FloatMatrix +calc_distances (float (*func)(short, short), const boolNDArray& bw, + short *xdist, short *ydist) { FloatMatrix dist (bw.dims ()); edtfunc (func, bw, xdist, ydist); const int numel = dist.numel (); + float* dist_vec = dist.fortran_vec (); for (int i = 0; i < numel; i++) - dist(i) = (*func)(xdist[i], ydist[i]); + dist_vec[i] = (*func)(xdist[i], ydist[i]); return dist; } template <class T> -T calc_index (Matrix bw, short *xdist, short *ydist) +T calc_index (const boolNDArray& bw, const short *xdist, const short *ydist) { + typedef typename T::element_type P; + T idx (bw.dims ()); const int numel = bw.numel (); const int rows = bw.rows (); + + P* idx_vec = idx.fortran_vec (); for(int i = 0; i < numel; i++) - idx (i) = i+1 - xdist[i] - ydist[i]*rows; + idx_vec[i] = i+1 - xdist[i] - ydist[i]*rows; + return idx; } @@ -472,10 +479,10 @@ // for matlab compatibility, we do not actually check if the values are all // 0 and 1, any non-zero value is considered true - const Matrix bw = args (0).matrix_value (); + const boolNDArray bw = args (0).bool_array_value (); if (error_state) { - error ("bwdist: BW must be a matrix"); + error ("bwdist: BW must be a logical matrix"); return retval; } @@ -512,9 +519,10 @@ dist = calc_distances (euclidean, bw, xdist, ydist); const Array<octave_idx_type> positions = (!bw).find (); const int zpos = positions.numel(); - for (int i = 0; i < zpos; i++) { - dist (positions(i)) = sqrt(dist(positions(i))); - } + const octave_idx_type* pos_vec = positions.fortran_vec (); + float* dist_vec = dist.fortran_vec (); + for (int i = 0; i < zpos; i++) + dist_vec[pos_vec[i]] = sqrt (dist_vec[pos_vec[i]]); } else if (method == "chessboard") dist = calc_distances (chessboard, bw, xdist, ydist); @@ -540,7 +548,7 @@ } /* -%!shared bw, out +%!shared bw %! %! bw = [0 1 0 1 0 1 1 0 %! 0 0 0 1 1 0 0 0 @@ -550,7 +558,8 @@ %! 1 1 1 1 0 0 0 1 %! 1 1 1 0 0 0 1 0 %! 0 0 1 0 0 0 1 1]; -%! + +%!test %! out = [ 1.00000 0.00000 1.00000 0.00000 1.00000 0.00000 0.00000 1.00000 %! 1.41421 1.00000 1.00000 0.00000 0.00000 1.00000 1.00000 1.41421 %! 2.23607 2.00000 1.00000 0.00000 0.00000 1.00000 2.00000 2.00000 @@ -561,10 +570,11 @@ %! 1.00000 1.00000 0.00000 1.00000 2.00000 1.00000 0.00000 0.00000]; %! out = single (out); %! -%!assert (bwdist (bw), out, 0.0001); # default is euclidean -%!assert (bwdist (bw, "euclidean"), out, 0.0001); -%!assert (bwdist (logical (bw), "euclidean"), out, 0.0001); -%! +%! assert (bwdist (bw), out, 0.0001); # default is euclidean +%! assert (bwdist (bw, "euclidean"), out, 0.0001); +%! assert (bwdist (logical (bw), "euclidean"), out, 0.0001); + +%!test %! out = [ 1 0 1 0 1 0 0 1 %! 1 1 1 0 0 1 1 1 %! 2 2 1 0 0 1 2 2 @@ -575,8 +585,9 @@ %! 1 1 0 1 2 1 0 0]; %! out = single (out); %! -%!assert (bwdist (bw, "chessboard"), out); -%! +%! assert (bwdist (bw, "chessboard"), out); + +%!test %! out = [ 1 0 1 0 1 0 0 1 %! 2 1 1 0 0 1 1 2 %! 3 2 1 0 0 1 2 2 @@ -587,8 +598,9 @@ %! 1 1 0 1 2 1 0 0]; %! out = single (out); %! -%!assert (bwdist (bw, "cityblock"), out); -%! +%! assert (bwdist (bw, "cityblock"), out); + +%!test %! out = [ 1.00000 0.00000 1.00000 0.00000 1.00000 0.00000 0.00000 1.00000 %! 1.41421 1.00000 1.00000 0.00000 0.00000 1.00000 1.00000 1.41421 %! 2.41421 2.00000 1.00000 0.00000 0.00000 1.00000 2.00000 2.00000 @@ -599,12 +611,35 @@ %! 1.00000 1.00000 0.00000 1.00000 2.00000 1.00000 0.00000 0.00000]; %! out = single (out); %! -%!assert (bwdist (bw, "quasi-euclidean"), out, 0.0001); +%! assert (bwdist (bw, "quasi-euclidean"), out, 0.0001); %! %! bw(logical (bw)) = 3; # there is no actual check if matrix is binary or 0 and 1 -%!assert (bwdist (bw, "quasi-euclidean"), out, 0.0001); +%! assert (bwdist (bw, "quasi-euclidean"), out, 0.0001); +%! %! bw(logical (bw)) = -2; # anything non-zero is considered object -%!assert (bwdist (bw, "quasi-euclidean"), out, 0.0001); +%! assert (bwdist (bw, "quasi-euclidean"), out, 0.0001); + +%!test +%! bw = [ 1 1 1 1 0 1 1 1 1 +%! 1 1 1 1 0 1 1 1 1 +%! 1 1 0 1 1 1 1 1 1 +%! 0 1 1 1 1 1 1 1 1]; %! -%!error bwdist (bw, "not a valid method"); +%! dist = [ 0 0 0 0 1 0 0 0 0 +%! 0 0 0 0 1 0 0 0 0 +%! 0 0 1 0 0 0 0 0 0 +%! 1 0 0 0 0 0 0 0 0]; +%! dist = single (dist); +%! +%! c = [ 1 5 9 13 13 21 25 29 33 +%! 2 6 10 14 14 22 26 30 34 +%! 3 7 10 15 19 23 27 31 35 +%! 8 8 12 16 20 24 28 32 36]; +%! c = uint32 (c); +%! +%! [dout, cout] = bwdist (bw, "euclidean"); +%! assert (dout, dist) +%! assert (cout, c) + +%!error <unknown METHOD> bwdist (bw, "not a valid method"); */