From: <ha...@us...> - 2008-03-06 11:48:20
|
Revision: 4698 http://octave.svn.sourceforge.net/octave/?rev=4698&view=rev Author: hauberg Date: 2008-03-06 03:48:17 -0800 (Thu, 06 Mar 2008) Log Message: ----------- Added radon transformation Modified Paths: -------------- trunk/octave-forge/main/image/INDEX Added Paths: ----------- trunk/octave-forge/main/image/inst/radon.m Modified: trunk/octave-forge/main/image/INDEX =================================================================== --- trunk/octave-forge/main/image/INDEX 2008-03-04 11:20:47 UTC (rev 4697) +++ trunk/octave-forge/main/image/INDEX 2008-03-06 11:48:17 UTC (rev 4698) @@ -47,6 +47,7 @@ stretchlim makelut applylut deriche + radon Black and white image functions bwarea bwdist Added: trunk/octave-forge/main/image/inst/radon.m =================================================================== --- trunk/octave-forge/main/image/inst/radon.m (rev 0) +++ trunk/octave-forge/main/image/inst/radon.m 2008-03-06 11:48:17 UTC (rev 4698) @@ -0,0 +1,88 @@ +## Copyright (C) 2007 Alexander Barth <bar...@gm...> +## +## +## 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 Foundation; either version 3 of the License, or (at +## your option) any later version. +## +## This program is distributed in the hope that it will be useful, but +## WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +## General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with this program; see the file COPYING. If not, see +## <http://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## @deftypefn {Function File} {[@var{RT},@var{xp}] =} radon(@var{I}, @var{theta}) +## @deftypefnx {Function File} {[@var{RT},@var{xp}] =} radon(@var{I}) +## +## Calculates the 2D-Radon transform of the matrix @var{I} at angles given in +## @var{theta}. To each element of @var{theta} corresponds a column in @var{RT}. +## The variable @var{xp} represents the x-axis of the rotated coordinate. +## If @var{theta} is not defined, then 0:179 is assumed. +## @end deftypefn + + +function [RT,xp] = radon (I,theta) + + ## Input checking + if (nargin == 0 | nargin > 2) + print_usage (); + elseif (nargin == 1) + theta = 0:179; + endif + + if (!ismatrix(I) || ndims(I) != 2) + error("radon: first input must be a MxN matrix"); + endif + + if (!isvector(theta)) + error("radon: second input must be a vector"); + endif + + [m, n] = size (I); + + # center of image + xc = floor ((m+1)/2); + yc = floor ((n+1)/2); + + # divide each pixel into 2x2 subpixels + + d = reshape (I,[1 m 1 n]); + d = d([1 1],:,[1 1],:); + d = reshape (d,[2*m 2*n])/4; + + b = ceil (sqrt (sum (size (I).^2))/2 + 1); + xp = [-b:b]'; + sz = size(xp); + + [X,Y] = ndgrid (0.75 - xc + [0:2*m-1]/2,0.75 - yc + [0:2*n-1]/2); + + X = X(:)'; + Y = Y(:)'; + d = d(:)'; + + th = theta*pi/180; + + for l=1:length (theta) + # project each pixel to vector (-sin(th),cos(th)) + Xp = -sin (th(l)) * X + cos (th(l)) * Y; + + ip = Xp + b + 1; + + k = floor (ip); + frac = ip-k; + + RT(:,l) = accumarray (k',d .* (1-frac),sz) + accumarray (k'+1,d .* frac,sz); + endfor + +endfunction + +%!test +%! A = radon(ones(2,2),30); +%! assert (A,[0 0 0.608253175473055 2.103325780167649 1.236538105676658 0.051882938682637 0]',1e-10) + + This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <ha...@us...> - 2008-03-11 08:52:01
|
Revision: 4708 http://octave.svn.sourceforge.net/octave/?rev=4708&view=rev Author: hauberg Date: 2008-03-11 01:52:02 -0700 (Tue, 11 Mar 2008) Log Message: ----------- Rewrite of 'edge'. This includes the addition of the 'canny' edge detector Modified Paths: -------------- trunk/octave-forge/main/image/INDEX trunk/octave-forge/main/image/inst/edge.m trunk/octave-forge/main/image/src/Makefile Added Paths: ----------- trunk/octave-forge/main/image/src/nonmax_supress.cc Modified: trunk/octave-forge/main/image/INDEX =================================================================== --- trunk/octave-forge/main/image/INDEX 2008-03-09 21:29:01 UTC (rev 4707) +++ trunk/octave-forge/main/image/INDEX 2008-03-11 08:52:02 UTC (rev 4708) @@ -48,6 +48,7 @@ makelut applylut deriche radon + nonmax_supress Black and white image functions bwarea bwdist Modified: trunk/octave-forge/main/image/inst/edge.m =================================================================== --- trunk/octave-forge/main/image/inst/edge.m 2008-03-09 21:29:01 UTC (rev 4707) +++ trunk/octave-forge/main/image/inst/edge.m 2008-03-11 08:52:02 UTC (rev 4708) @@ -1,60 +1,128 @@ +## Copyright (C) 2008 S\xF8ren Hauberg +## +## 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 Foundation; either version 3 of the License, or +## (at your option) any later version. +## +## This program is distributed in the hope that it will be useful, +## but WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +## GNU General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with this program; if not, write to the Free Software +## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +## +## Note: The implementation and help text for the 'andy' edge detector came with +## the following notice: +## Copyright (C) 1999 Andy Adler +## This code has no warrany whatsoever. +## Do what you like with this code as long as you +## leave this copyright in place. + ## -*- texinfo -*- -## @deftypefn {Function File} {[@var{imout}, @var{thresh}] = } edge(@var{im}, @var{method}, @var{thresh}, @var{param2}) -## Find image edges. +## @deftypefn {Function File} @var{bw} = edge (@var{im}) +## @deftypefnx {Function File} @var{bw} = edge (@var{im}, "sobel") +## Finds the edges in @var{im} using the Sobel approximation to the +## derivatives. Edge points are defined as points where the length of +## the gradient exceeds a threshold and is larger than it's neighbours +## in either the horizontal or vertical direction. +## @deftypefnx {Function File} @var{bw} = edge (@var{im}, "sobel", @var{thresh}) +## Same as above except @var{thresh} is used as a threshold. +## @deftypefnx {Function File} @var{bw} = edge (@var{im}, "sobel", @var{thresh}, @var{direction}) +## Same as above except the derivate only is approximated in @var{direction}, +## where @var{direction} can be either "horizontal", "vertical", or "both" +## (default). +## @deftypefnx {Function File} [@var{bw}, @var{thresh} ] = edge (@var{im}, "sobel", ...) +## Same as any of the above except the used threshold is alose returned. ## -## The output is -## @table @code -## @item @var{imout} -## Output image -## @item @var{thresh} -## Output thresholds -## @end table +## @deftypefnx {Function File} @var{bw} = edge (@var{im}, "prewitt") +## Finds the edges in @var{im} using the Prewitt approximation to the +## derivatives. Edge points are defined as points where the length of +## the gradient exceeds a threshold and is larger than it's neighbours +## in either the horizontal or vertical direction. +## @deftypefnx {Function File} @var{bw} = edge (@var{im}, "prewitt", @var{thresh}) +## Same as above except @var{thresh} is used as a threshold. +## @deftypefnx {Function File} @var{bw} = edge (@var{im}, "prewitt", @var{thresh}, @var{direction}) +## Same as above except the derivate only is approximated in @var{direction}, +## where @var{direction} can be either "horizontal", "vertical", or "both" +## (default). +## @deftypefnx {Function File} [@var{bw}, @var{thresh} ] = edge (@var{im}, "prewitt", ...) +## Same as any of the above except the used threshold is alose returned. ## -## The input is -## @table @code -## @item @var{im} -## Input image (greyscale) -## @item @var{thresh} -## Threshold value (value is estimated if not given) -## @end table +## @deftypefnx {Function File} @var{bw} = edge (@var{im}, "roberts") +## Finds the edges in @var{im} using the Roberts approximation to the +## derivatives. Edge points are defined as points where the length of +## the gradient exceeds a threshold and is larger than it's neighbours +## in either the horizontal or vertical direction. +## @deftypefnx {Function File} @var{bw} = edge (@var{im}, "roberts", @var{thresh}) +## Same as above except @var{thresh} is used as a threshold. +## @deftypefnx {Function File} @var{bw} = edge (@var{im}, "roberts", @var{thresh}, @var{option}) +## Same as above except the thinning step can be turned off by setting +## @var{option} to "nothinning". The default value for @var{option} is +## "thinning". +## @deftypefnx {Function File} [@var{bw}, @var{thresh} ] = edge (@var{im}, "roberts", ...) +## Same as any of the above except the used threshold is alose returned. +## +## @deftypefnx {Function File} @var{bw} = edge (@var{im}, "kirsch") +## Finds the edges in @var{im} using the Kirsch approximation to the +## derivatives. Edge points are defined as points where the length of +## the gradient exceeds a threshold and is larger than it's neighbours +## in either the horizontal or vertical direction. +## @deftypefnx {Function File} @var{bw} = edge (@var{im}, "kirsch", @var{thresh}) +## Same as above except @var{thresh} is used as a threshold. +## @deftypefnx {Function File} @var{bw} = edge (@var{im}, "kirsch", @var{thresh}, @var{direction}) +## Same as above except the derivate only is approximated in @var{direction}, +## where @var{direction} can be either "horizontal", "vertical", or "both" +## (default). +## @deftypefnx {Function File} [@var{bw}, @var{thresh} ] = edge (@var{im}, "kirsch", ...) +## Same as any of the above except the used threshold is alose returned. ## -## The following methods are based on high pass filtering the image in -## two directions, calculating a combined edge weight from and then thresholding +## @deftypefnx{Function File} @var{bw} = edge(@var{im}, "log") +## Finds edges in @var{im} by convolving with the Laplacian of Gaussian +## filter, and finding zero crossings. Only zero crossings where the +## filter response is larger than an automaticly computed threshold. +## @deftypefnx{Function File} @var{bw} = edge(@var{im}, "log", @var{thresh}) +## Same as above except @var{thresh} is used as threshold instead of the +## automaticly computed threshold. +## @deftypefnx{Function File} @var{bw} = edge(@var{im}, "log", @var{thresh}, @var{sigma}) +## Same as above, except the spread of the Laplacian of Gaussian is given +## by @var{sigma}. The default value is 2. +## @deftypefnx{Function File} [@var{bw}, @var{thresh}] = edge(@var{im}, "log", ...) +## Same as any of the above, except the used threshold is also returned. +## +## @deftypefnx{Function File} @var{bw} = edge(@var{im}, "zerocross", @var{thresh}, @var{filter}) +## Finds edges in the image @var{im} by convolving it with the user-supplied filter +## @var{filter} and finding zero crossings larger than @var{thresh}. If @var{thresh} is +## [] a threshold will be chosen automaticly. +## @deftypefnx{Function File} [@var{bw}, @var{thresh}] = edge(@var{im}, "zerocross", ...) +## Same as above, except the used threshold is returned. ## -## @table @code -## @item method = 'roberts' -## @example -## filt1= [1 0 ; 0 -1]; filt2= rot90( filt1 ) -## combine= sqrt( filt1^2 + filt2^2 ) -## @end example -## @item method = 'sobel' -## @example -## filt1= [1 2 1;0 0 0;-1 -2 -1]; filt2= rot90( filt1 ) -## combine= sqrt( filt1^2 + filt2^2 ) -## @end example -## @item method = 'prewitt' -## @example -## filt1= [1 1 1;0 0 0;-1 -1 -1]; filt2= rot90( filt1 ) -## combine= sqrt( filt1^2 + filt2^2 ) -## @end example -## @item method = 'kirsh' -## @example -## filt1= [1 2 1;0 0 0;-1 -2 -1]; filt2 .. filt8 are 45 degree rotations of filt1 -## combine= max( filt1 ... filt8 ) -## @end example -## @end table +## @deftypefnx{Function File} @var{bw} = edge(@var{im}, "canny") +## Finds edges using the Canny edge detector. +## @deftypefnx{Function File} @var{bw} = edge(@var{im}, "canny", @var{thresh}) +## Finds edges using the Canny edge detector, where the thresholds used in the +## hysteresis are user-defined. If @var{thresh} is a two dimensional vector it's first +## element is used as the lower threshold, while the second element is used as the +## high threshold. If, on the other hand, @var{thresh} is a single scalar it is used as +## the high threshold, while the lower threshold is 0.4*@var{thresh}. +## @deftypefnx{Function File} @var{bw} = edge(@var{im}, "canny", @var{thresh}, @var{sigma}) +## Same as above except the spread of the low-pass filtering Gaussian is given by +## @var{sigma} (defaults to 2). +## @deftypefnx{Function File} [@var{bw}, @var{threshold}] = edge(@var{im}, "canny", ...) +## Same as any of the above except except the two used thresholds are returned as a +## vector in @var{threshold}. ## -## Methods based on filtering the image and finding zero crossings +## @deftypefnx{Function File} @var{bw} = edge(@var{im}, "lindeberg") +## Finds edges using in @var{im} using the differential geomtric single-scale edge +## detector given by Tony Lindeberg. +## @deftypefnx{Function File} @var{bw} = edge(@var{im}, "lindeberg", @var{sigma}) +## Same as above except the scale can be controlled with the parameter @var{sigma} +## (defaults to 2). ## -## @table @code -## @item method = 'log' -## Laplacian of Gaussian. -## @var{param2} is the standard deviation of the filter, default is 2. -## @item method = 'zerocross' -## generic zero-crossing filter. -## @var{param2} is the user supplied filter. -## @item method = 'andy' -## A.Adler's idea (c) 1999. somewhat based on the canny method. The steps are +## @deftypefnx{Function File} @var{bw} = edge(@var{im}, "andy") +## A.Adler's idea (c) 1999. Somewhat based on the canny method. The steps are ## @enumerate ## @item ## Do a sobel edge detection and to generate an image at @@ -70,125 +138,365 @@ ## the HT image. ## @end enumerate ## -## The parameters for the method are: -## @table @code -## @item param2(1)==0 or 4 or 8 +## The parameters for the method is given in a vector: +## @table @asis +## @item params(1)==0 or 4 or 8 ## Perform x connected dilatation (step 3). -## @item param2(2) +## @item params(2) ## Dilatation coeficient (threshold) in step 3. -## @item param2(3) +## @item params(3) ## Length of edge extention convolution (step 2). -## @item param2(4) +## @item params(4) ## Coeficient of extention convolution in step 2. ## @end table -## defaults = [8 1 3 3] -## @end table +## defaults = [8, 1, 3, 3] +## +## @seealso{fspecial, nonmax_supress} ## @end deftypefn -# Copyright (C) 1999 Andy Adler -# This code has no warrany whatsoever. -# Do what you like with this code as long as you -# leave this copyright in place. -# -# $Id$ +function [bw, out_threshold, g45_out, g135_out] = edge (im, method, varargin) + ## Get the image + if (nargin == 0) + error("edge: not enough input arguments"); + endif + if ( !isgray(im) ) + error("edge: first input must be a gray-scale image"); + endif -function [imout, thresh] = edge( im, method, thresh, param2 ) + ## Get the method + if (nargin == 1) + method = "Sobel"; + endif + if (!ischar(method)) + error("edge: second argument must be a string"); + endif + method = lower(method); -[n,m]= size(im); -xx= 2:m-1; -yy= 2:n-1; + ## Perform the actual edge detection + switch (method) + ##################################### + ## S O B E L + ##################################### + case "sobel" + ## Get the direction argument + direction = get_direction(varargin{:}); + ## Create filters; + h1 = fspecial("sobel"); # horizontal + h3 = h1'; # vertical + ## Compute edge strength + switch(direction) + case "horizontal" + strength = abs( conv2(im, h1, "same") ); + case "vertical" + strength = abs( conv2(im, h3, "same") ); + case "both" + strength = sqrt( conv2(im, h1, "same").^2 + ... + conv2(im, h3, "same").^2 ); + endswitch + ## Get threshold + if (nargin > 2 && isscalar(varargin{1})) + thresh = varargin{1}; + else + thresh = mean(abs(strength(:))); + endif + ## Perform thresholding and simple thinning + strength(strength<=thresh) = 0; + bw = simple_thinning(strength); -if strcmp(method,'roberts') || strcmp(method,'sobel') || ... - strcmp(method,'prewitt') - + ##################################### + ## P R E W I T T + ##################################### + case "prewitt" + ## Get the direction argument + direction = get_direction(varargin{:}); + ## Create filters; + h1 = fspecial("prewitt"); # vertical + h3 = h1'; # horizontal + ## Compute edge strength + switch(direction) + case "vertical" + strength = abs( conv2(im, h1, "same") ); + case "horizontal" + strength = abs( conv2(im, h3, "same") ); + case "both" + strength = sqrt( conv2(im, h1, "same").^2 + ... + conv2(im, h3, "same").^2 ); + endswitch + ## Get threshold + if (nargin > 2 && isscalar(varargin{1})) + thresh = varargin{1}; + else + thresh = mean(abs(strength(:))); + endif + ## Perform thresholding and simple thinning + strength(strength<=thresh) = 0; + bw = simple_thinning(strength); + + ##################################### + ## K I R S C H + ##################################### + case "kirsch" + ## Get the direction argument + direction = get_direction(varargin{:}); + ## Create filters; + h1 = fspecial("kirsch"); # vertical + h3 = h1'; # horizontal + ## Compute edge strength + switch(direction) + case "vertical" + strength = abs( conv2(im, h1, "same") ); + case "horizontal" + strength = abs( conv2(im, h3, "same") ); + case "both" + strength = sqrt( conv2(im, h1, "same").^2 + ... + conv2(im, h3, "same").^2 ); + endswitch + ## Get threshold + if nargin > 2 && isscalar(varargin{1}) + thresh = varargin{1}; + else + thresh = mean(abs(strength(:))); + endif + ## Perform thresholding and simple thinning + strength(strength<=thresh) = 0; + bw = simple_thinning(strength); - if strcmp(method,'roberts') - filt= [1 0;0 -1]/4; tv= 6; - elseif strcmp(method,'sobel') - filt= [1 2 1;0 0 0; -1 -2 -1]/8; tv= 2; - elseif strcmp(method,'prewitt') - filt= [1 1 1;0 0 0; -1 -1 -1]/6; tv= 4; - end + ##################################### + ## R O B E R T S + ##################################### + case "roberts" + ## Get the thinning argument (option) + if (nargin == 4) + option = varargin{2}; + if (!ischar(option)) + error("edge: 'option' must be a string"); + endif + option = lower(option); + if (!any(strcmp(option, {"thinning", "nothinning"}))) + error("edge: 'option' must be either 'thinning', or 'nothinning'"); + endif + else + option = "thinning"; + endif + ## Create filters; + h1 = [1 0; 0 -1]; + h2 = [0 1; -1 0]; + ## Compute edge strength + g45 = conv2(im, h1, "same"); + g135 = conv2(im, h2, "same"); + strength = abs( g45 ) + abs( g135 ); + ## Get threshold + if (nargin > 2 && isscalar(varargin{1})) + thresh = varargin{1}; + else + thresh = mean(abs(strength(:))); + endif + ## Perform thresholding and simple thinning + strength(strength<=thresh) = 0; + if strcmp(option, "thinning") + bw = simple_thinning(strength); + else + bw = (strength > 0); + endif + ## Check if g45 and g135 should be returned + if (nargout == 4) + g45_out = g45; + g135_out = g135; + endif + + ##################################### + ## L A P L A C I A N O F G A U S S I A N + ##################################### + case "log" + ## Get sigma + if (nargin == 4 && isscalar(varargin{2})) + sigma = varargin{2}; + else + sigma = 2; + endif + ## Create the filter + s = ceil(3*sigma); + %[x y] = meshgrid(-s:s); + %f = (x.^2 + y.^2 - sigma^2) .* exp(-(x.^2 + y.^2)/(2*sigma^2)); + %f = f/sum(f(:)); + f = fspecial("log", 2*s+1, sigma); + ## Perform convolution with the filter f + g = conv2(im, f, "same"); + ## Get threshold + if (nargin > 2 && isscalar(varargin{1})) + thresh = varargin{1}; + else + thresh = mean(abs(g(:))); + endif + ## Find zero crossings + zc = zerocrossings(g); + bw = (abs(g) >= thresh) & zc; + + ##################################### + ## Z E R O C R O S S I N G + ##################################### + case "zerocross" + ## Get the filter + if (nargin == 4 && ismatrix(varargin{2})) + f = varargin{2}; + else + error("edge: a filter must be given as the fourth argument when 'zerocross' is used"); + endif + ## Perform convolution with the filter f + g = conv2(im, f, "same"); + ## Get threshold + if (nargin > 2 && isscalar(varargin{1})) + thresh = varargin{1}; + else + thresh = mean(abs(g(:))); + endif + ## Find zero crossings + zc = zerocrossings(g); + bw = (abs(g) >= thresh) & zc; - imo= conv2(im, rot90(filt), 'same').^2 + conv2(im, filt, 'same').^2; - -# check to see if the user supplied a threshold -# if not, calculate one in the same way as Matlab + ##################################### + ## C A N N Y + ##################################### + case "canny" + ## Get sigma + if (nargin == 4 && isscalar(varargin{2})) + sigma = varargin{2}; + else + sigma = 2; + endif + ## Change scale + gauss = fspecial("gaussian", 2*ceil(3*sigma)+1, sigma); + J = conv2(im, gauss, "same"); + ## Canny enhancer + p = [1 0 -1]/2; + Jx = conv2(J, p, "same"); + Jy = conv2(J, p', "same"); + Es = sqrt( Jx.^2 + Jy.^2 ); + Eo = atan2(Jy,Jx); + ## Get thresholds + if (nargin > 2 && isscalar(varargin{1})) + thresh = [0.4*varargin{1}, varargin{1}]; + elseif (nargin > 2 && ismatrix(varargin{1}) && len(varargin{1}(:)) == 2) + thresh = varargin{1}(:); + else + tmp = mean(abs(Es(:))); + thresh = [0.4*tmp, tmp]; + endif + bw = nonmax_supress(Es, Eo, thresh(1), thresh(2)); - if nargin<3 - thresh= sqrt( tv* mean(mean( imo(yy,xx) )) ); - end + ##################################### + ## L I N D E B E R G + ##################################### + case "lindeberg" + ## In case the user asks for more then 1 output argument + ## we define thresh to be -1. + thresh = -1; + ## Get sigma + if (nargin > 2 && isscalar(varargin{1})) + sigma = varargin{1}; + else + sigma = 2; + endif + ## Filters for computing the derivatives + Px = [-1 0 1; -1 0 1; -1 0 1]; + Py = [1 1 1; 0 0 0; -1 -1 -1]; + Pxx = conv2(Px, Px, "full"); + Pyy = conv2(Py, Py, "full"); + Pxy = conv2(Px, Py, "full"); + Pxxx = conv2(Pxx, Px, "full"); + Pyyy = conv2(Pyy, Py, "full"); + Pxxy = conv2(Pxx, Py, "full"); + Pxyy = conv2(Pyy, Px, "full"); + ## Change scale + gauss = fspecial("gaussian", 2*ceil(3*sigma)+1, sigma); + L = conv2(im, gauss, "same"); + ## Compute derivatives + Lx = conv2(L, Px, "same"); + Ly = conv2(L, Py, "same"); + Lxx = conv2(L, Pxx, "same"); + Lyy = conv2(L, Pyy, "same"); + Lxy = conv2(L, Pxy, "same"); + Lxxx = conv2(L, Pxxx, "same"); + Lyyy = conv2(L, Pyyy, "same"); + Lxxy = conv2(L, Pxxy, "same"); + Lxyy = conv2(L, Pxyy, "same"); + ## Compute directional derivatives + Lvv = Lx.^2.*Lxx + 2.*Lx.*Ly.*Lxy + Ly.^2.*Lyy; + Lvvv = Lx.^3.*Lxxx + 3.*Lx.^2.*Ly.*Lxxy ... + + 3.*Lx.*Ly.^2.*Lxyy + 3.*Ly.^3.*Lyyy; + ## Perform edge detection + bw = zerocrossings(Lvv) & Lvvv < 0; -# The filters are defined for sqrt(imo), but since we calculated imo, compare -# to thresh ^2 + ##################################### + ## A N D Y + ##################################### + case "andy" + [bw, out_threshold] = andy (im, method, varargin{:}); + + otherwise + error("edge: unsupported edge detector: %s", method); + endswitch + + if (nargout > 1) + out_threshold = thresh; + endif +endfunction - imout= ( imo >= thresh^2 ); +## An auxilary function that parses the 'direction' argument from 'varargin' +function direction = get_direction(varargin) + if (nargin >= 2) + direction = varargin{2}; + if (!ischar(direction)) + error("edge: direction must be a string"); + endif + direction = lower(direction); + if (!any(strcmp(direction, {"horizontal", "vertical", "both"}))) + error("edge :direction must be either 'horizontal', 'vertical', or 'both'"); + endif + else + direction = "both"; + endif +endfunction -# Thin the wide edges - xpeak= imo(yy,xx-1) <= imo(yy,xx) & imo(yy,xx) > imo(yy,xx+1) ; - ypeak= imo(yy-1,xx) <= imo(yy,xx) & imo(yy,xx) > imo(yy+1,xx) ; - imout(yy,xx)= imout(yy,xx) & ( xpeak | ypeak ); +## An auxilary function that performs a very simple thinning. +## Strength is an image containing the edge strength. +## bw contains a 1 in (r,c) if +## 1) strength(r,c) is greater than both neighbours in the +## vertical direction, OR +## 2) strength(r,c) is greater than both neighbours in the +## horizontal direction. +## Note the use of OR. +function bw = simple_thinning(strength) + [r c] = size(strength); + x = ( strength > [ zeros(r,1) strength(:,1:end-1) ] & ... + strength > [ strength(:,2:end) zeros(r,1) ] ); + y = ( strength > [ zeros(1,c); strength(1:end-1,:) ] & ... + strength > [ strength(2:end,:); zeros(1,c) ] ); + bw = x | y; +endfunction -elseif strcmp(method,'kirsch') +## Auxilary function. Finds the zero crossings of the +## 2-dimensional function f. (By Etienne Grossmann) +function z = zerocrossings(f) + z0 = f<0; ## Negative + [R,C] = size(f); + z = zeros(R,C); + z(1:R-1,:) |= z0(2:R,:); ## Grow + z(2:R,:) |= z0(1:R-1,:); + z(:,1:C-1) |= z0(:,2:C); + z(:,2:C) |= z0(:,1:C-1); - filt1= [1 2 1;0 0 0;-1 -2 -1]; fim1= conv2(im,filt1,'same'); - filt2= [2 1 0;1 0 -1;0 -1 -2]; fim2= conv2(im,filt2,'same'); - filt3= [1 0 -1;2 0 -2;1 0 -1]; fim3= conv2(im,filt3,'same'); - filt4= [0 1 2;-1 0 1;-2 -1 0]; fim4= conv2(im,filt4,'same'); + z &= !z0; ## "Positive zero-crossings"? +endfunction - imo= reshape(max([abs(fim1(:)) abs(fim2(:)) abs(fim3(:)) abs(fim4(:))]'),n,m); +## The 'andy' edge detector that was present in older versions of 'edge'. +## The function body has simply been copied from the old implementation. +## -- Soren Hauberg, march 11th, 2008 +function [imout, thresh] = andy(im, method, thresh, param2) + [n,m]= size(im); + xx= 2:m-1; + yy= 2:n-1; - if nargin<3 - thresh= 2* mean(mean( imo(yy,xx) )) ; - end - - imout= imo >= thresh ; - -# Thin the wide edges - xpeak= imo(yy,xx-1) <= imo(yy,xx) & imo(yy,xx) > imo(yy,xx+1) ; - ypeak= imo(yy-1,xx) <= imo(yy,xx) & imo(yy,xx) > imo(yy+1,xx) ; - imout(yy,xx)= imout(yy,xx) & ( xpeak | ypeak ); - -elseif strcmp(method,'log') || strcmp(method,'zerocross') - - if strcmp(method,'log') - if nargin >= 4; sd= param2; - else sd= 2; - end - - sz= ceil(sd*3); - [x,y]= meshgrid( -sz:sz, -sz:sz ); - filt = exp( -( x.^2 + y.^2 )/2/sd^2 ) .* ... - ( x.^2 + y.^2 - 2*sd^2 ) / 2 / pi / sd^6 ; - else - filt = param2; - end - filt = filt - mean(filt(:)); - - imo= conv2(im, filt, 'same'); - - if nargin<3 || isempty( thresh ) - thresh= 0.75* mean(mean( abs(imo(yy,xx)) )) ; - end - - zcross= imo > 0; - yd_zc= diff( zcross ); - xd_zc= diff( zcross' )'; - yd_io= abs(diff( imo ) ) > thresh; - xd_io= abs(diff( imo')') > thresh; - -# doing it this way puts the transition at the <=0 point - xl= zeros(1,m); yl= zeros(n,1); - imout= [ ( yd_zc == 1 ) & yd_io ; xl] | ... - [xl; ( yd_zc == -1 ) & yd_io ] | ... - [ ( xd_zc == 1 ) & xd_io , yl] | ... - [yl, ( xd_zc == -1 ) & xd_io ]; - -elseif strcmp(method,'canny') - error("method canny not implemented"); - -elseif strcmp(method,'andy') - filt= [1 2 1;0 0 0; -1 -2 -1]/8; tv= 2; imo= conv2(im, rot90(filt), 'same').^2 + conv2(im, filt, 'same').^2; if nargin<3 || thresh==[]; @@ -244,42 +552,4 @@ % imout = bwselect( imee, rem(ff-1, n)+1, ceil(ff/n), 8); imout = imee; - -else - - error (['Method ' method ' is not recognized']); - -end - - - - -# -# $Log$ -# Revision 1.3 2007/01/04 23:41:47 hauberg -# Minor changes in help text -# -# Revision 1.2 2007/01/02 21:58:38 hauberg -# Documentation is now in Texinfo (looks better on the website) -# -# Revision 1.1 2006/08/20 12:59:32 hauberg -# Changed the structure to match the package system -# -# Revision 1.1 2002/03/17 02:38:52 aadler -# fill and edge detection operators -# -# Revision 1.4 2000/11/20 17:13:07 aadler -# works? -# -# Revision 1.3 1999/06/09 17:29:36 aadler -# implemented 'andy' mode edge detection -# -# Revision 1.2 1999/06/08 14:26:50 aadler -# zero-cross and LoG filters work -# -# Revision 1.1 1999/06/07 21:01:38 aadler -# Initial revision -# -# - -# +endfunction Modified: trunk/octave-forge/main/image/src/Makefile =================================================================== --- trunk/octave-forge/main/image/src/Makefile 2008-03-09 21:29:01 UTC (rev 4707) +++ trunk/octave-forge/main/image/src/Makefile 2008-03-11 08:52:02 UTC (rev 4708) @@ -13,7 +13,7 @@ endif all: cordflt2.oct bwlabel.oct bwfill.oct rotate_scale.oct \ - houghtf.oct graycomatrix.oct deriche.oct __bwdist.oct \ + houghtf.oct graycomatrix.oct deriche.oct __bwdist.oct nonmax_supress.oct \ $(JPEG) $(PNG) $(IMAGEMAGICK) jpgread.oct: jpgread.cc Added: trunk/octave-forge/main/image/src/nonmax_supress.cc =================================================================== --- trunk/octave-forge/main/image/src/nonmax_supress.cc (rev 0) +++ trunk/octave-forge/main/image/src/nonmax_supress.cc 2008-03-11 08:52:02 UTC (rev 4708) @@ -0,0 +1,188 @@ +/* +Copyright (C) 2005 S\xF8ren Hauberg + +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 Foundation; either version 2 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with this program; if not, write to the Free Software +Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + +Author: S\xF8ren Hauberg <hauberg at gmail dot com> + +2005-04-16 S\xF8ren Hauberg <hauberg at gmail dot com> +* Initial revision + +*/ + +#include <octave/oct.h> +#include <math.h> +#include <stack> +#include <utility> + +static bool any_bad_argument(const octave_value_list& args) { + /* TODO: We should check that the matrices have the same size. */ + int nargin = args.length(); + if ( (nargin >= 2) && args(0).is_real_matrix() && args(1).is_real_matrix() ) + { + if (nargin == 2) + { + return false; + } else if (nargin == 4 && args(2).is_real_scalar() && + args(3).is_real_scalar() ) + { + return false; + } + } + error("Input arguments must be two 2-dimensional matrices of the same size."); + return true; +} + +DEFUN_DLD(nonmax_supress,args,nargout,"\ +-*- texinfo -*-\n\ +@deftypefn {Function File} nonmax_supress (@var{Es}, @var{Eo})\n\ +Performs non-maximum supression on the given edge data. \ +@var{Es} is a matrix containing the edge strength (the length of \ +the gradient), and @var{Eo} is the edge normal orientation (the \ +direction of the gradient).\n\ +\n\ +@deftypefnx {Function File} nonmax_supress (@var{Es}, @var{Eo},\ + @var{low}, @var{high} )\n\ +Performs non-maximum supression and hysteresis thresholdong, using \ +@var{low} and @var{high} as thresholds.\n\ +\n\ +This function is designed to be used as part of the Canny edge \ +detection, and not to be used in general. So if you use this function: \ +Beware...\n\ +\n\ +@seealso{edge}\n\ +@end deftypefn\n\ +") +{ + octave_value_list retval; + if (any_bad_argument(args)) { + return retval; + } + std::stack< std::pair<int,int> > S; + + /* Neighbourhood directions in radians */ + const double d[4] = { + 0.0, + M_PI * 45.0 / 180.0, + M_PI * 90.0 / 180.0, + M_PI * 135.0 / 180.0 + }; + + const Matrix Es = args(0).matrix_value(); + Matrix Eo = args(1).matrix_value(); + double low, high; + bool hysteresis = (args.length()==4); + if (hysteresis) { + low = args(2).scalar_value(); + high = args(3).scalar_value(); + } else { + low = high = 0; + } + + const int rows = Es.rows(); + const int cols = Es.columns(); + + /**************************** + ** Non-maximum supression ** + ****************************/ + Matrix In = Matrix( rows, cols, 0.0 ); + for (int r = 1; r < rows-1; r++) { + for (int c = 1; c < cols-1; c++) { + const double orientation = Eo(r,c); + const double strength = Es(r,c); + + int best_d = 0; + const double dist = fabs( orientation-d[0] ); + for (int i = 1; i < 4; i++) { + if ( fabs( orientation-d[i] ) < dist ) { best_d = i; } + } + Eo(r,c) = best_d; + + switch (best_d) { + case 0: // 0 degrees + if ( (strength > Es(r,c-1)) && (strength > Es(r,c+1)) ) + { In(r,c) = strength; } + break; + case 1: // 45 degrees + if ( (strength > Es(r-1,c+1)) && (strength > Es(r+1,c-1)) ) + { In(r,c) = strength; } + break; + case 2: // 90 degrees + if ( (strength > Es(r-1,c)) && (strength > Es(r+1,c)) ) + { In(r,c) = strength; } + break; + case 3: // 135 degrees + if ( (strength > Es(r-1,c-1)) && (strength > Es(r+1,c+1)) ) + { In(r,c) = strength; } + break; + } + + if (hysteresis && In(r,c) > high) { + S.push( std::pair<int,int>(r,c) ); + } + } + } + + if (hysteresis == false) { + retval.append(In); + return retval; + } + + /************************** + ** Hysteresis threshold ** + **************************/ + boolMatrix out = boolMatrix( rows, cols, false ); + while (S.empty() == false) { + std::pair<int, int> p = S.top(); + S.pop(); + const int r = p.first; + const int c = p.second; + if (r < 0 || r >= rows || c < 0 || c >= cols || out(r,c) == true) + { continue; } + + out(r,c) = true; + const int dir = (int)Eo(r,c); + switch (dir) { + case 0: // 0 degrees + if ( In(r-1,c) > low ) + { S.push(std::pair<int,int>(r-1,c)); } + if ( In(r+1,c) > low ) + { S.push(std::pair<int,int>(r+1,c)); } + break; + case 1: // 45 degrees + if ( In(r-1,c-1) > low ) + { S.push(std::pair<int,int>(r-1,c-1)); } + if ( In(r+1,c+1) > low ) + { S.push(std::pair<int,int>(r+1,c+1)); } + break; + case 2: // 90 degrees + if ( In(r,c-1) > low ) + { S.push(std::pair<int,int>(r,c-1)); } + if ( In(r,c+1) > low ) + { S.push(std::pair<int,int>(r,c+1)); } + break; + case 3: // 135 degrees + if ( In(r-1,c+1) > low ) + { S.push(std::pair<int,int>(r-1,c+1)); } + if ( In(r+1,c-1) > low ) + { S.push(std::pair<int,int>(r+1,c-1)); } + break; + } + } + + retval.append(out); + return retval; +} + This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <ha...@us...> - 2008-03-21 14:01:33
|
Revision: 4770 http://octave.svn.sourceforge.net/octave/?rev=4770&view=rev Author: hauberg Date: 2008-03-21 07:01:36 -0700 (Fri, 21 Mar 2008) Log Message: ----------- New function for locating spatial maximas Modified Paths: -------------- trunk/octave-forge/main/image/INDEX Added Paths: ----------- trunk/octave-forge/main/image/inst/immaximas.m Modified: trunk/octave-forge/main/image/INDEX =================================================================== --- trunk/octave-forge/main/image/INDEX 2008-03-21 12:19:52 UTC (rev 4769) +++ trunk/octave-forge/main/image/INDEX 2008-03-21 14:01:36 UTC (rev 4770) @@ -33,6 +33,7 @@ graycomatrix houghtf graythresh + immaximas Filtering imfilter colfilt Added: trunk/octave-forge/main/image/inst/immaximas.m =================================================================== --- trunk/octave-forge/main/image/inst/immaximas.m (rev 0) +++ trunk/octave-forge/main/image/inst/immaximas.m 2008-03-21 14:01:36 UTC (rev 4770) @@ -0,0 +1,101 @@ +## Copyright (c) 2003-2005 Peter Kovesi +## School of Computer Science & Software Engineering +## The University of Western Australia +## http://www.csse.uwa.edu.au/ +## +## Permission is hereby granted, free of charge, to any person obtaining a copy +## of this software and associated documentation files (the "Software"), to deal +## in the Software without restriction, subject to the following conditions: +## +## The above copyright notice and this permission notice shall be included in all +## copies or substantial portions of the Software. +## +## The Software is provided "as is", without warranty of any kind. +## +## I've made minor changes compared to the original 'nonmaxsuppts' function developed +## by Peter Kovesi. The original is available at +## http://www.csse.uwa.edu.au/~pk/research/matlabfns/Spatial/nonmaxsuppts.m +## -- Søren Hauberg, 2008 + +## -*- texinfo -*- +## @deftypefn {Function File} {[@var{r}, @var{c}] =} immaximas (@var{im}, @var{radius}) +## @deftypefnx{Function File} {[@var{r}, @var{c}] =} immaximas (@var{im}, @var{radius}, @var{thresh}) +## Finds local spatial maximas of the given image. A local spatial maxima is +## defined as an image point with a value that is larger than all neighbouring +## values in a square region of width 2*@var{radius}+1. By default @var{radius} +## is 1, such that a 3 by 3 neighbourhood is searched. If the @var{thresh} input +## argument is supplied, only local maximas with a value greater than @var{thresh} +## are retained. +## +## The output vectors @var{r} and @var{c} contain the row-column coordinates +## of the local maximas. The actual values are computed to sub-pixel precision +## by fitting a parabola to the data around the pixel. +## @end deftypefn + +function [r, c] = immaximas(im, radius, thresh) + ## Check input + if (nargin == 0) + error("immaximas: not enough input arguments"); + endif + if (nargin <= 1 || isempty(radius)) + radius = 1; + endif + if (nargin <= 2) + thresh = []; + endif + if (!ismatrix(im) || ndims(im) != 2) + error("immaximas: first input argument must be an M by N matrix"); + endif + if (!isscalar(radius) && !isempty(radius)) + error("immaximas: second input argument must be a scalar or an empty matrix"); + endif + if (!isscalar(thresh) && !isempty(thresh)) + error("immaximas: third input argument must be a scalar or an empty matrix"); + endif + + ## Find local maximas + s = size(im); + sze = 2*radius+1; + mx = ordfilt2(im, sze^2, ones(sze, "logical")); + mx2 = ordfilt2(im, sze^2-1, ones(sze, "logical")); + + ## Make mask to exclude points within radius of the image boundary. + bordermask = zeros(s, "logical"); + bordermask(radius+1:end-radius, radius+1:end-radius) = 1; + + # Find maxima, threshold, and apply bordermask + immx = (im == mx) & (im != mx) & bordermask; + if (!isempty(thresh)) + immx &= (im>thresh); + endif + + ## Find local maximas and fit parabolas locally + [r, c] = find(immx); + if (!isempty(r)) + ind = sub2ind(s,r,c); # 1D indices of feature points + w = 1; # Width that we look out on each side of the feature point to fit a local parabola + + ## Indices of points above, below, left and right of feature point + indrminus1 = max(ind-w,1); + indrplus1 = min(ind+w,s(1)*s(2)); + indcminus1 = max(ind-w*s(1),1); + indcplus1 = min(ind+w*s(1),s(1)*s(2)); + + ## Solve for quadratic down rows + cy = im(ind); + ay = (im(indrminus1) + im(indrplus1))/2 - cy; + by = ay + cy - im(indrminus1); + rowshift = -w*by./(2*ay); # Maxima of quadradic + + ## Solve for quadratic across columns + cx = im(ind); + ax = (im(indcminus1) + im(indcplus1))/2 - cx; + bx = ax + cx - im(indcminus1); + colshift = -w*bx./(2*ax); # Maxima of quadradic + + ## Add subpixel corrections to original row and column coords. + r += rowshift; + c += colshift; + endif +endfunction + This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <ha...@us...> - 2008-03-21 15:42:11
|
Revision: 4776 http://octave.svn.sourceforge.net/octave/?rev=4776&view=rev Author: hauberg Date: 2008-03-21 08:41:49 -0700 (Fri, 21 Mar 2008) Log Message: ----------- Added N-dimensional version of 'ordfilt2' Modified Paths: -------------- trunk/octave-forge/main/image/INDEX Added Paths: ----------- trunk/octave-forge/main/image/inst/ordfiltn.m Modified: trunk/octave-forge/main/image/INDEX =================================================================== --- trunk/octave-forge/main/image/INDEX 2008-03-21 15:11:39 UTC (rev 4775) +++ trunk/octave-forge/main/image/INDEX 2008-03-21 15:41:49 UTC (rev 4776) @@ -43,7 +43,8 @@ imadjust imnoise medfilt2 - ordfilt2 cordflt2 + ordfilt2 + ordfiltn uintlut stretchlim makelut applylut Added: trunk/octave-forge/main/image/inst/ordfiltn.m =================================================================== --- trunk/octave-forge/main/image/inst/ordfiltn.m (rev 0) +++ trunk/octave-forge/main/image/inst/ordfiltn.m 2008-03-21 15:41:49 UTC (rev 4776) @@ -0,0 +1,99 @@ +## Copyright (C) 2008 Soren Hauberg +## +## 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 Foundation; either version 3 +## of the License, or (at your option) any later version. +## +## This program is distributed in the hope that it will be useful, but +## WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +## General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with this program; If not, see <http://www.gnu.org/licenses/>. + +## This function is based on 'ordfiltn' by Teemu Ikonen which is released under +## GPLv2 or later. + +## -*- texinfo -*- +## @deftypefn {Function File} {} ordfiltn(@var{A}, @var{nth}, @var{domain}, [@var{S}, @var{padding}]) +## Two dimensional ordered filtering. +## +## Ordered filter replaces an element of @var{A} with the @var{nth} +## element of the sorted set of neighbours defined by the logical +## (boolean) matrix @var{domain}. +## Neighbour elements are selected to the sort if the corresponding +## element in the @var{domain} matrix is true. +## +## The optional variable @var{S} is a matrix of size(@var{domain}). +## Values of @var{S} corresponding to nonzero values of domain are +## added to values obtained from @var{A} when doing the sorting. +## +## Optional variable @var{padding} determines how the matrix @var{A} +## is padded from the edges. See @code{padarray} for details. +## +## @seealso{ordfilt2, padarray} +## @end deftypefn + + +## Author: Teemu Ikonen <tpi...@pc...> +## Created: 5.5.2000 +## Keywords: image processing filtering + +function retval = ordfiltn(A, nth, domain, varargin) + ## Check input + if (nargin < 3) + error("ordfiltn: not enough input arguments"); + endif + if (!ismatrix(A)) + error("ordfiltn: first input must be an array"); + endif + if (!isscalar(nth) || nth <= 0 || nth != round(nth)) + error("ordfiltn: second input argument must be a positive integer"); + endif + if (!ismatrix(domain) && !isscalar(domain)) + error("ordfiltn: third input argument must be an array or a scalar"); + endif + if (isscalar(domain) && (domain <= 0 || domain != round(domain))) + error("ordfiltn: third input argument must be a positive integer, when it is a scalar"); + endif + if (isscalar(domain)) + domain = ones(repmat(domain, 1, ndims(A)), "logical"); + endif + + if (ndims(A) != ndims(domain)) + error("ordfiltn: first and second argument must have same dimensionality"); + endif + if (any(size(A) < size(domain))) + error("ordfiltn: domain array cannot be larger than the data array"); + endif + + ## Parse varargin + S = zeros(size(domain)); + padding = 0; + for i=1:length(varargin) + a = varargin{:}; + if (ischar(a) || isscalar(a)) + padding = a; + elseif (ismatrix(a) && size_equal(a, domain)) + S = a; + endif + endfor + + ## Make sure 'domain' is logical. The C++ code assumes this. + domain = logical(domain); + + ## Pad array + pad = floor((size(domain)+1)/2); + A = padarray(A, pad, padding); + even = ( round(size(domain)/2) == size(domain)/2 ); + idx = cell(1, ndims(A)); + for k = 1:ndims(A) + idx{k} = (even(k)+1):size(A,k); + endfor + A = A(idx{:}); + + ## Perform the filtering + retval = __cordfltn__(A, nth, domain, S); +endfunction This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <ha...@us...> - 2008-03-23 23:14:40
|
Revision: 4786 http://octave.svn.sourceforge.net/octave/?rev=4786&view=rev Author: hauberg Date: 2008-03-23 16:14:44 -0700 (Sun, 23 Mar 2008) Log Message: ----------- Removed 'cordflt2.cc' as '__cordfltn__' is more general. Modified Paths: -------------- trunk/octave-forge/main/image/src/Makefile Added Paths: ----------- trunk/octave-forge/main/image/inst/cordflt2.m Removed Paths: ------------- trunk/octave-forge/main/image/src/cordflt2.cc Added: trunk/octave-forge/main/image/inst/cordflt2.m =================================================================== --- trunk/octave-forge/main/image/inst/cordflt2.m (rev 0) +++ trunk/octave-forge/main/image/inst/cordflt2.m 2008-03-23 23:14:44 UTC (rev 4786) @@ -0,0 +1,27 @@ +## Copyright (C) 2008 Soren Hauberg +## +## 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 Foundation; either version 3 of the License, or +## (at your option) any later version. +## +## This program is distributed in the hope that it will be useful, +## but WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +## GNU General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with this program; If not, see <http://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## @deftypefn {Function File} cordflt2(@var{A}, @var{nth}, @var{domain}, @var{S}) +## Implementation of two-dimensional ordered filtering. This function has been +## deprecated and should NOT be used. Instead use @code{ordfilt2}. +## @seealso{ordfilt2} +## @end deftypefn + +function varargout = cordflt2(varargin) + warning(["cordflt2: this function is deprecated and will be removed in upcoming " + "releases. Use 'ordfilt2' instead."]); + [varargout{1:nargout}] = __cordfltn__(varargin{:}); +endfunction Modified: trunk/octave-forge/main/image/src/Makefile =================================================================== --- trunk/octave-forge/main/image/src/Makefile 2008-03-23 21:29:37 UTC (rev 4785) +++ trunk/octave-forge/main/image/src/Makefile 2008-03-23 23:14:44 UTC (rev 4786) @@ -12,7 +12,7 @@ IMAGEMAGICK=__magick_read__.oct endif -all: cordflt2.oct __cordfltn__.oct bwlabel.oct bwfill.oct rotate_scale.oct \ +all: __cordfltn__.oct bwlabel.oct bwfill.oct rotate_scale.oct \ hough_line.oct graycomatrix.oct deriche.oct __bwdist.oct nonmax_supress.oct \ $(JPEG) $(PNG) $(IMAGEMAGICK) Deleted: trunk/octave-forge/main/image/src/cordflt2.cc =================================================================== --- trunk/octave-forge/main/image/src/cordflt2.cc 2008-03-23 21:29:37 UTC (rev 4785) +++ trunk/octave-forge/main/image/src/cordflt2.cc 2008-03-23 23:14:44 UTC (rev 4786) @@ -1,211 +0,0 @@ -// Copyright (C) 2000 Teemu Ikonen -// -// 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 Foundation; either version 2 -// of the License, or (at your option) any later version. -// -// This program is distributed in the hope that it will be useful, but -// WITHOUT ANY WARRANTY; without even the implied warranty of -// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU -// General Public License for more details. -// -// You should have received a copy of the GNU General Public License -// along with this program; if not, see <http://www.gnu.org/licenses/>. - -#include <octave/oct.h> - -#define SWAP(a, b) { SWAP_temp = (a); (a)=(b); (b) = SWAP_temp; } - -// Template function for comparison -// ET is the type of the matrix element -template <class ET> -inline bool compare(const ET a, const ET b) -{ - if(a > b) - return 1; - else - return 0; -} - -// Explicit template function for complex compare -template <> inline bool compare<Complex>(const Complex a, const Complex b) -{ - double anorm2 = a.real() * a.real() + a.imag() * a.imag(); - double bnorm2 = b.real() * b.real() + b.imag() * b.imag(); - - if( anorm2 > bnorm2 ) { - return 1; - } else { - return 0; - } -} - -// select nth largest member from the array values -// Partitioning algorithm, see Numerical recipes chap. 8.5 -template <class ET> -ET selnth(ET *vals, int len, int nth) -{ - ET SWAP_temp; - ET hinge; - int l, r, mid, i, j; - - l = 0; - r = len - 1; - for(;;) { - // if partition size is 1 or two, then sort and return - if(r <= l+1) { - if(r == l+1 && compare<ET>(vals[l], vals[r])) { - SWAP(vals[l], vals[r]); - } - return vals[nth]; - } else { - mid = (l+r) >> 1; - SWAP(vals[mid], vals[l+1]); - // choose median of l, mid, r to be the hinge element - // and set up sentinels in the borders (order l, l+1 and r) - if(compare<ET>(vals[l], vals[r])) { - SWAP(vals[l], vals[r]); - } - if(compare<ET>(vals[l+1], vals[r])) { - SWAP(vals[l+1], vals[r]); - } - if(compare<ET>(vals[l], vals[l+1])) { - SWAP(vals[l], vals[l+1]); - } - i = l+1; - j = r; - hinge = vals[l+1]; - for(;;) { - do i++; while(compare<ET>(hinge, vals[i])); - do j--; while(compare<ET>(vals[j], hinge)); - if(i > j) - break; - SWAP(vals[i], vals[j]); - } - vals[l+1] = vals[j]; - vals[j] = hinge; - if(j >= nth) - r = j - 1; - if(j <= nth) - l = i; - } - } -} - -// Template function for doing the actual filtering -// MT is the type of the matrix to be filtered (Matrix or ComplexMatrix) -// ET is the type of the element of the matrix (double or Complex) -template <class MT, class ET> -octave_value_list do_filtering(MT A, int nth, boolMatrix dom, MT S) -{ - int i, j, c, d; - - int len = 0; - for(j = 0; j < dom.columns(); j++) { - for(i = 0; i < dom.rows(); i++) { - if(dom.elem(i,j)) - len++; - } - } - if(nth > len - 1) { - warning("nth should be less than number of non-zero values in domain"); - warning("setting nth to largest possible value\n"); - nth = len - 1; - } - if(nth < 0) { - warning("nth should be non-negative, setting to 1\n"); - nth = 0; // nth is a c-index - } - - int rowoffset = (dom.columns() + 1)/2 - 1; - int coloffset = (dom.rows() + 1)/2 - 1; - - //outputs - octave_value_list out; - const int origx = A.columns() - dom.columns()+1; - const int origy = A.rows() - dom.rows()+1; - MT retval = MT(origy, origx); - - int *offsets = new int[len]; - ET *values = new ET[len]; - ET *adds = new ET[len]; - - c = 0; - d = A.rows(); - for(j = 0; j < dom.columns(); j++) { - for(i = 0; i < dom.rows(); i++) { - if(dom.elem(i,j)) { - offsets[c] = (i - coloffset) + (j - rowoffset)*d; - adds[c] = S.elem(i,j); - c++; - } - } - } - - ET *data = A.fortran_vec(); - int base = coloffset + A.rows()*rowoffset; - for(j = 0; j < retval.columns(); j++) { - for(i = 0; i < retval.rows(); i++) { - for(c = 0; c < len; c++) { - values[c] = data[base + offsets[c]] + adds[c]; - } - base++; - retval(i, j) = selnth(values, len, nth); - } - base += dom.rows() - 1; - } - - out(0) = octave_value(retval); - - return out; -} - -// instantiate template functions -template bool compare<double>(const double, const double); -template double selnth(double *, int, int); -template Complex selnth(Complex *, int, int); -template octave_value_list do_filtering<Matrix, double>(Matrix, int, boolMatrix, Matrix); -// g++ is broken, explicit instantiation of specialized template function -// confuses the compiler. -//template int compare<Complex>(const Complex, const Complex); -template octave_value_list do_filtering<ComplexMatrix, Complex>(ComplexMatrix, int, boolMatrix, ComplexMatrix); - -DEFUN_DLD(cordflt2, args, , "\ --*- texinfo -*-\n\ -@deftypefn {Function File} cordflt2(@var{A}, @var{nth}, @var{domain}, @var{S})\n\ -Implementation of two-dimensional ordered filtering. In general this function\n\ -should NOT be used. Instead use @code{ordfilt2}.\n\ -@seealso{ordfilt2}\n\ -@end deftypefn\n\ -") -{ - if(args.length() != 4) { - print_usage (); - return octave_value_list(); - } - - // nth is an index to an array, thus - 1 - int nth = (int) (args(1).vector_value())(0) - 1; - boolMatrix dom = args(2).bool_matrix_value(); - - octave_value_list retval; - - if(args(0).is_real_matrix()) { - Matrix A = args(0).matrix_value(); - Matrix S = args(3).matrix_value(); - retval = do_filtering<Matrix, double>(A, nth, dom, S); - } - else if(args(0).is_complex_matrix()) { - ComplexMatrix A = args(0).complex_matrix_value(); - ComplexMatrix S = args(3).complex_matrix_value(); - retval = do_filtering<ComplexMatrix, Complex>(A, nth, dom, S); - } - else { - error("A should be real or complex matrix\n"); - return octave_value_list(); - } - - return retval; - -} This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <ha...@us...> - 2008-03-24 11:25:56
|
Revision: 4788 http://octave.svn.sourceforge.net/octave/?rev=4788&view=rev Author: hauberg Date: 2008-03-24 04:26:00 -0700 (Mon, 24 Mar 2008) Log Message: ----------- More morphological operators Modified Paths: -------------- trunk/octave-forge/main/image/INDEX Added Paths: ----------- trunk/octave-forge/main/image/inst/bwhitmiss.m trunk/octave-forge/main/image/inst/imclose.m trunk/octave-forge/main/image/inst/imdilate.m trunk/octave-forge/main/image/inst/imerode.m trunk/octave-forge/main/image/inst/imopen.m trunk/octave-forge/main/image/inst/imtophat.m Modified: trunk/octave-forge/main/image/INDEX =================================================================== --- trunk/octave-forge/main/image/INDEX 2008-03-24 10:52:45 UTC (rev 4787) +++ trunk/octave-forge/main/image/INDEX 2008-03-24 11:26:00 UTC (rev 4788) @@ -64,6 +64,13 @@ erode bwborder edge conndef + bwhitmiss +Morhophological Operations + imerode + imdilate + imopen + imclose + imtophat Colour controls cmpermute cmunique Added: trunk/octave-forge/main/image/inst/bwhitmiss.m =================================================================== --- trunk/octave-forge/main/image/inst/bwhitmiss.m (rev 0) +++ trunk/octave-forge/main/image/inst/bwhitmiss.m 2008-03-24 11:26:00 UTC (rev 4788) @@ -0,0 +1,64 @@ +## Copyright (C) 2008 Soren Hauberg +## +## 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 Foundation; either version 2 +## of the License, or (at your option) any later version. +## +## This program is distributed in the hope that it will be useful, +## but WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +## GNU General Public License for more details. + +## -*- texinfo -*- +## @deftypefn {Function File} @var{bw2} = bwhitmiss (@var{bw1}, @var{se1}, @var{se1}) +## @deftypefnx{Function File} @var{bw2} = bwhitmiss (@var{bw1}, @var{interval}) +## Perform the binary hit-miss operation. +## +## If two structuring elements @var{se1} and @var{se1} are given, the hit-miss +## operation is defined as +## @example +## bw2 = erode(bw1, se1) & erode(!bw1, se2); +## @end example +## If instead an 'interval' array is given, two structuring elements are computed +## as +## @example +## se1 = (interval == 1) +## se2 = (interval == -1) +## @end example +## and then the operation is defined as previously. +## @seealso{bwmorph} +## @end deftypefn + +function bw = bwhitmiss(im, varargin) + ## Checkinput + if (nargin != 2 && nargin != 3) + print_usage(); + endif + if (!ismatrix(im) || !isreal(im)) + error("bwhitmiss: first input argument must be a real matrix"); + endif + + ## Get structuring elements + if (nargin == 2) # bwhitmiss (im, interval) + interval = varargin{1}; + if (!isreal(interval)) + error("bwhitmiss: second input argument must be a real matrix"); + endif + if (!all( (interval(:) == 1) | (interval(:) == 0) | (interval(:) == -1) )) + error("bwhitmiss: second input argument can only contain the values -1, 0, and 1"); + endif + se1 = (interval == 1); + se2 = (interval == -1); + else # bwhitmiss (im, se1, se2) + se1 = varargin{1}; + se2 = varargin{2}; + if (!all((se1(:) == 1) | (se1(:) == 0)) || !all((se2(:) == 1) | (se2(:) == 0))) + error("bwhitmiss: structuring elements can only contain zeros and ones."); + endif + endif + + ## Perform filtering + bw = erode(im, se1) & erode(!im, se2); + +endfunction Added: trunk/octave-forge/main/image/inst/imclose.m =================================================================== --- trunk/octave-forge/main/image/inst/imclose.m (rev 0) +++ trunk/octave-forge/main/image/inst/imclose.m 2008-03-24 11:26:00 UTC (rev 4788) @@ -0,0 +1,41 @@ +## Copyright (C) 2008 Soren Hauberg +## +## 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 Foundation; either version 2 +## of the License, or (at your option) any later version. +## +## This program is distributed in the hope that it will be useful, +## but WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +## GNU General Public License for more details. + +## -*- texinfo -*- +## @deftypefn {Function File} @var{B} = imclose (@var{A}, @var{se}) +## Perform morphological closening on a given image. +## The image @var{A} must be a grayscale or binary image, and @var{se} must be a +## structuring element. +## +## The closening corresponds to a dilation followed by an erosion of the image, i.e. +## @example +## B = imerode(imdilate(A, se), se); +## @end example +## @seealso{imdilate, imerode, imclose} +## @end deftypefn + +function retval = imclose(im, se) + ## Checkinput + if (nargin != 2) + print_usage(); + endif + if (!ismatrix(im) || !isreal(im)) + error("imclose: first input argument must be a real matrix"); + endif + if (!ismatrix(se) || !isreal(se)) + error("imclose: second input argument must be a real matrix"); + endif + + ## Perform filtering + retval = imerode(imdilate(im, se), se); + +endfunction Added: trunk/octave-forge/main/image/inst/imdilate.m =================================================================== --- trunk/octave-forge/main/image/inst/imdilate.m (rev 0) +++ trunk/octave-forge/main/image/inst/imdilate.m 2008-03-24 11:26:00 UTC (rev 4788) @@ -0,0 +1,35 @@ +## Copyright (C) 2008 Soren Hauberg +## +## 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 Foundation; either version 3 +## of the License, or (at your option) any later version. +## +## This program is distributed in the hope that it will be useful, +## but WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +## GNU General Public License for more details. + +## -*- texinfo -*- +## @deftypefn {Function File} @var{B} = imdilate (@var{A}, @var{se}) +## Perform morphological dilation on a given image. +## The image @var{A} must be a grayscale or binary image, and @var{se} must be a +## structuring element. +## @seealso{imerode, imopen, imclose} +## @end deftypefn + +function retval = imdilate(im, se) + ## Checkinput + if (nargin != 2) + print_usage(); + endif + if (!ismatrix(im) || !isreal(im)) + error("imdilate: first input argument must be a real matrix"); + endif + if (!ismatrix(se) || !isreal(se)) + error("imdilate: second input argument must be a real matrix"); + endif + + ## Perform filtering + retval = ordfiltn(im, sum(se(:)!=0), se, 0); +endfunction Added: trunk/octave-forge/main/image/inst/imerode.m =================================================================== --- trunk/octave-forge/main/image/inst/imerode.m (rev 0) +++ trunk/octave-forge/main/image/inst/imerode.m 2008-03-24 11:26:00 UTC (rev 4788) @@ -0,0 +1,35 @@ +## Copyright (C) 2008 Soren Hauberg +## +## 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 Foundation; either version 3 +## of the License, or (at your option) any later version. +## +## This program is distributed in the hope that it will be useful, +## but WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +## GNU General Public License for more details. + +## -*- texinfo -*- +## @deftypefn {Function File} @var{B} = imerode (@var{A}, @var{se}) +## Perform morphological erosion on a given image. +## The image @var{A} must be a grayscale or binary image, and @var{se} must be a +## structuring element. +## @seealso{imdilate, imopen, imclose} +## @end deftypefn + +function retval = imerode(im, se) + ## Checkinput + if (nargin != 2) + print_usage(); + endif + if (!ismatrix(im) || !isreal(im)) + error("imerode: first input argument must be a real matrix"); + endif + if (!ismatrix(se) || !isreal(se)) + error("imerode: second input argument must be a real matrix"); + endif + + ## Perform filtering + retval = ordfiltn(im, 1, se, 0); +endfunction Added: trunk/octave-forge/main/image/inst/imopen.m =================================================================== --- trunk/octave-forge/main/image/inst/imopen.m (rev 0) +++ trunk/octave-forge/main/image/inst/imopen.m 2008-03-24 11:26:00 UTC (rev 4788) @@ -0,0 +1,41 @@ +## Copyright (C) 2008 Soren Hauberg +## +## 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 Foundation; either version 2 +## of the License, or (at your option) any later version. +## +## This program is distributed in the hope that it will be useful, +## but WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +## GNU General Public License for more details. + +## -*- texinfo -*- +## @deftypefn {Function File} @var{B} = imopen (@var{A}, @var{se}) +## Perform morphological opening on a given image. +## The image @var{A} must be a grayscale or binary image, and @var{se} must be a +## structuring element. +## +## The opening corresponds to an erosion followed by a dilation of the image, i.e. +## @example +## B = imdilate(imerode(A, se), se); +## @end example +## @seealso{imdilate, imerode, imclose} +## @end deftypefn + +function retval = imopen(im, se) + ## Checkinput + if (nargin != 2) + print_usage(); + endif + if (!ismatrix(im) || !isreal(im)) + error("imopen: first input argument must be a real matrix"); + endif + if (!ismatrix(se) || !isreal(se)) + error("imopen: second input argument must be a real matrix"); + endif + + ## Perform filtering + retval = imdilate(imerode(im, se), se); + +endfunction Added: trunk/octave-forge/main/image/inst/imtophat.m =================================================================== --- trunk/octave-forge/main/image/inst/imtophat.m (rev 0) +++ trunk/octave-forge/main/image/inst/imtophat.m 2008-03-24 11:26:00 UTC (rev 4788) @@ -0,0 +1,35 @@ +## Copyright (C) 2005 Carvalho-Mariel +## +## 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 Foundation; either version 2 +## of the License, or (at your option) any later version. +## +## This program is distributed in the hope that it will be useful, +## but WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +## GNU General Public License for more details. + +## -*- texinfo -*- +## @deftypefn {Function File} @var{B} = imtophat (@var{A}, @var{se}) +## Perform morphological top hat filtering. +## The image @var{A} must be a grayscale or binary image, and @var{se} must be a +## structuring element. +## @end deftypefn + +function retval = imtophat(im, se) + ## Checkinput + if (nargin != 2) + print_usage(); + endif + if (!ismatrix(im) || !isreal(im)) + error("imtophat: first input argument must be a real matrix"); + endif + if (!ismatrix(se) || !isreal(se)) + error("imtophat: second input argument must be a real matrix"); + endif + + ## Perform filtering + retval = im & !imopen(im, SE); + +endfunction This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <ha...@us...> - 2008-03-24 11:26:56
|
Revision: 4789 http://octave.svn.sourceforge.net/octave/?rev=4789&view=rev Author: hauberg Date: 2008-03-24 04:26:58 -0700 (Mon, 24 Mar 2008) Log Message: ----------- added rgbplot Modified Paths: -------------- trunk/octave-forge/main/image/INDEX Added Paths: ----------- trunk/octave-forge/main/image/inst/rgbplot.m Modified: trunk/octave-forge/main/image/INDEX =================================================================== --- trunk/octave-forge/main/image/INDEX 2008-03-24 11:26:00 UTC (rev 4788) +++ trunk/octave-forge/main/image/INDEX 2008-03-24 11:26:58 UTC (rev 4789) @@ -3,6 +3,7 @@ image imagesc imshow + rgbplot Read/write imread imwrite Added: trunk/octave-forge/main/image/inst/rgbplot.m =================================================================== --- trunk/octave-forge/main/image/inst/rgbplot.m (rev 0) +++ trunk/octave-forge/main/image/inst/rgbplot.m 2008-03-24 11:26:58 UTC (rev 4789) @@ -0,0 +1,34 @@ +## Copyright (C) 2005 Berge-Gladel +## +## 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 Foundation; either version 2 +## of the License, or (at your option) any later version. +## +## This program is distributed in the hope that it will be useful, +## but WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +## GNU General Public License for more details. + +## -*- texinfo -*- +## @deftypefn {Function File} rgbplot (@var{map}) +## @deftypefnx{Function File} @var{h} = rgbplot (@var{map}) +## Plot a given color map. +## The matrix @var{map} must be a @math{M} by 3 matrix. The three columns of the +## colormap matrix are plotted in red, green, and blue lines. +## +## If an output is requested, a graphics handle to the plot is returned. +## @end deftypefn + +function h_out = rgbplot(map) + ## Check input + if (!ismatrix(map) || ndims(map) != 2 || columns(map) != 3) + error("rgbplot: input must be a M by 3 matrix"); + endif + + ## Plot + h = plot(map(i,1), "-r", map(i,2), "g-", map(i,3), "b-"); + if (nargout > 0) + h_out = h; + endif +endfunction This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <ha...@us...> - 2008-03-24 11:57:29
|
Revision: 4790 http://octave.svn.sourceforge.net/octave/?rev=4790&view=rev Author: hauberg Date: 2008-03-24 04:57:21 -0700 (Mon, 24 Mar 2008) Log Message: ----------- new function: imcomplement Modified Paths: -------------- trunk/octave-forge/main/image/INDEX Added Paths: ----------- trunk/octave-forge/main/image/inst/imcomplement.m Modified: trunk/octave-forge/main/image/INDEX =================================================================== --- trunk/octave-forge/main/image/INDEX 2008-03-24 11:26:58 UTC (rev 4789) +++ trunk/octave-forge/main/image/INDEX 2008-03-24 11:57:21 UTC (rev 4790) @@ -96,6 +96,7 @@ rgb2gray rgb2ind label2rgb + imcomplement Colour maps flag lines Added: trunk/octave-forge/main/image/inst/imcomplement.m =================================================================== --- trunk/octave-forge/main/image/inst/imcomplement.m (rev 0) +++ trunk/octave-forge/main/image/inst/imcomplement.m 2008-03-24 11:57:21 UTC (rev 4790) @@ -0,0 +1,45 @@ +## Copyright (C) 2008 Soren Hauberg +## +## 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 Foundation; either version 3 +## of the License, or (at your option) any later version. +## +## This program is distributed in the hope that it will be useful, but +## WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +## General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with 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. +## +## 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}}. +## @end deftypefn + +function B = imcomplement(A) + ## Check input + if (nargin != 1) + error("imcomplement: not enough input arguments"); + endif + if (!ismatrix(A)) + error("imcomplement: input must be an array"); + endif + + ## Take action depending on the class of A + if (isa(A, "double") || isa(A, "single")) + B = 1 - A; + elseif (islogical(A)) + B = !A; + elseif (isinteger(A)) + B = intmax(class(A)) - A; + else + error("imcomplement: unsupported input class: '%s'", class(A)); + endif +endfunction This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <ha...@us...> - 2008-03-24 19:16:35
|
Revision: 4799 http://octave.svn.sourceforge.net/octave/?rev=4799&view=rev Author: hauberg Date: 2008-03-24 12:15:21 -0700 (Mon, 24 Mar 2008) Log Message: ----------- imsmooth: support for bilateral filtering Modified Paths: -------------- trunk/octave-forge/main/image/inst/imsmooth.m trunk/octave-forge/main/image/src/Makefile Added Paths: ----------- trunk/octave-forge/main/image/src/__bilateral__.cc Modified: trunk/octave-forge/main/image/inst/imsmooth.m =================================================================== --- trunk/octave-forge/main/image/inst/imsmooth.m 2008-03-24 14:48:44 UTC (rev 4798) +++ trunk/octave-forge/main/image/inst/imsmooth.m 2008-03-24 19:15:21 UTC (rev 4799) @@ -31,6 +31,8 @@ ## Smoothing using a circular averaging linear filter. ## @item "Median" ## Median filtering. +## @item "Bilateral" +## Gaussian bilateral filtering. ## @item "Perona & Malik" ## @itemx "Perona and Malik" ## @itemx "P&M" @@ -76,6 +78,31 @@ ## ## The image is extrapolated symmetrically before the filtering is performed. ## +## @strong{Gaussian bilateral filtering} +## +## The image is smoothed using Gaussian bilateral filtering as described by +## Tomasi and Manduchi [2]. The filtering result is computed as +## @example +## @var{J(x0, y0) = k * SUM SUM @var{I}(x,y) * w(x, y, x0, y0, @var{I}(x0,y0), @var{I}(x,y)) +## x y +## @end example +## where @code{k} a normalisation variable, and +## @example +## w(x, y, x0, y0, @var{I}(x0,y0), @var{I}(x,y)) +## = exp(-0.5*d([x0,y0],[x,y])^2/@var{sigma_d}^2) +## * exp(-0.5*d(@var{I}(x0,y0),@var{I}(x,y))^2/@var{sigma_r}^2), +## @end example +## with @code{d} being the Euclidian distance function. The two paramteres +## @var{sigma_d} and @var{sigma_r} control the amount of smoothing. @var{sigma_d} +## is the size of the spatial smoothing filter, while @var{sigma_r} is the size +## of the range filter. When @var{sigma_r} is large the filter behaves almost +## like the isotropic Gaussian filter with spread @var{sigma_d}, and when it is +## small edges are preserved better. By default @var{sigma_d} is 2, and @var{sigma_r} +## is @math{10/255} for floating points images (with integer images this is +## multiplied with the maximal possible value representable by the integer class). +## +## The image is extrapolated symmetrically before the filtering is performed. +## ## @strong{Perona and Malik} ## ## The image is smoothed using anisotropic diffusion as described by Perona and @@ -126,6 +153,10 @@ ## IEEE Transactions on Pattern Analysis and Machine Intelligence, ## 12(7):629-639, 1990. ## +## [2] C. Tomasi and R. Manduchi, +## "Bilateral Filtering for Gray and Color Images", +## Proceedings of the 1998 IEEE International Conference on Computer Vision, Bombay, India. +## ## @seealso{imfilter, fspecial} ## @end deftypefn @@ -150,7 +181,6 @@ ## Save information for later C = class(I); - I = double(I); ## Take action depending on 'name' switch (lower(name)) @@ -171,7 +201,7 @@ h = ceil(3*s); f = exp( (-(-h:h).^2)./(2*s^2) ); f /= sum(f); ## Pad image - I = impad(I, h, h, "symmetric"); + I = double(impad(I, h, h, "symmetric")); ## Perform the filtering for i = imchannels:-1:1 J(:,:,i) = conv2(f, f, I(:,:,i), "valid"); @@ -196,7 +226,7 @@ f2 = ones(1,s(1))/s(1); f1 = ones(1,s(2))/s(2); ## Pad image - I = impad(I, floor([s(2), s(2)-1]/2), floor([s(1), s(1)-1]/2), "symmetric"); + I = impad(double(I), floor([s(2), s(2)-1]/2), floor([s(1), s(1)-1]/2), "symmetric"); ## Perform the filtering for i = imchannels:-1:1 J(:,:,i) = conv2(f1, f2, I(:,:,i), "valid"); @@ -218,7 +248,7 @@ ## Compute filter f = fspecial("disk", r); ## Pad image - I = impad(I, r, r, "symmetric"); + I = impad(double(I), r, r, "symmetric"); ## Perform the filtering for i = imchannels:-1:1 J(:,:,i) = conv2(I(:,:,i), f, "valid"); @@ -246,6 +276,36 @@ J(:,:,i) = medfilt2(I(:,:,i), s, "symmetric"); endfor + ############################### + ### Bilateral Filtering ### + ############################### + case "bilateral" + ## Check input + if (len > 0 && !isempty(varargin{1})) + if (isscalar(varargin{1}) && varargin{1} > 0) + sigma_d = varargin{1}; + else + error("imsmooth: spread of closeness function must be a positive scalar"); + endif + else + sigma_d = 2; + endif + if (len > 1 && !isempty(varargin{2})) + if (isscalar(varargin{2}) && varargin{2} > 0) + sigma_r = varargin{2}; + else + error("imsmooth: spread of similarity function must be a positive scalar"); + endif + else + sigma_r = 10/255; + if (isinteger(I)), sigma_r *= intmax(C); endif + endif + ## Pad image + s = max([round(3*sigma_d),1]); + I = impad(I, s, s, "symmetric"); + ## Perform the filtering + J = __bilateral__(I, sigma_d, sigma_r); + ############################ ### Perona and Malik ### ############################ @@ -291,6 +351,7 @@ endif endif ## Perform the filtering + I = double(I); for i = imchannels:-1:1 J(:,:,i) = pm(I(:,:,i), iter, lambda, method); endfor Modified: trunk/octave-forge/main/image/src/Makefile =================================================================== --- trunk/octave-forge/main/image/src/Makefile 2008-03-24 14:48:44 UTC (rev 4798) +++ trunk/octave-forge/main/image/src/Makefile 2008-03-24 19:15:21 UTC (rev 4799) @@ -12,7 +12,7 @@ IMAGEMAGICK=__magick_read__.oct endif -all: __cordfltn__.oct bwlabel.oct bwfill.oct rotate_scale.oct \ +all: __cordfltn__.oct __bilateral__.oct bwlabel.oct bwfill.oct rotate_scale.oct \ hough_line.oct graycomatrix.oct deriche.oct __bwdist.oct nonmax_supress.oct \ $(JPEG) $(PNG) $(IMAGEMAGICK) Added: trunk/octave-forge/main/image/src/__bilateral__.cc =================================================================== --- trunk/octave-forge/main/image/src/__bilateral__.cc (rev 0) +++ trunk/octave-forge/main/image/src/__bilateral__.cc 2008-03-24 19:15:21 UTC (rev 4799) @@ -0,0 +1,191 @@ +// Copyright (C) 2008 Soren Hauberg +// +// 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 Foundation; either version 3 +// of the License, or (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, but +// WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program; if not, see <http://www.gnu.org/licenses/>. + +#include <octave/oct.h> + +inline +int MAX(const int a, const int b) +{ + if (a > b) + return a; + else + return b; +} + +inline +double gauss1(const double x, const double mu, const double sigma) +{ + const double d = x-mu; + return exp( -0.5*(d*d)/(sigma*sigma) ); +} + +inline +double gauss(const double *x, const double *mu, const double sigma, const octave_idx_type ndims) +{ + double s = 0; + for (octave_idx_type i = 0; i < ndims; i++) + { + const double d = x[i] - mu[i]; + s += d*d; + } + return exp( -0.5*s/(sigma*sigma) ); +} + +template <class MatrixType> +octave_value +bilateral(const MatrixType &im, const double sigma_d, const double sigma_r) +{ + // Get sizes + const octave_idx_type ndims = im.ndims(); + const dim_vector size = im.dims(); + const octave_idx_type num_planes = (ndims == 2) ? 1 : size(2); + + // Build spatial kernel + const int s = MAX(round(3*sigma_d), 1); + Matrix kernel(2*s+1, 2*s+1); + for (octave_idx_type r = 0; r < 2*s+1; r++) + { + for (octave_idx_type c = 0; c < 2*s+1; c++) + { + const int dr = r-s; + const int dc = c-s; + kernel(r,c) = exp( -0.5*( dr*dr + dc*dc )/(sigma_d*sigma_d) ); + } + } + + // Allocate output + dim_vector out_size(size); + out_size(0) = MAX(size(0) - 2*s, 0); + out_size(1) = MAX(size(1) - 2*s, 0); + MatrixType out = MatrixType(out_size); + + // Iterate over every element of 'out'. + for (octave_idx_type r = 0; r < out_size(0); r++) + { + for (octave_idx_type c = 0; c < out_size(1); c++) + { + // For each neighbour + double val[num_planes]; + double sum[num_planes]; + double k = 0; + for (octave_idx_type i = 0; i < num_planes; i++) + { + val[i] = im(r,c,i); + sum[i] = 0; + } + for (octave_idx_type kr = 0; kr < 2*s+1; kr++) + { + for (octave_idx_type kc = 0; kc < 2*s+1; kc++) + { + double lval[num_planes]; + for (octave_idx_type i = 0; i < num_planes; i++) lval[i] = im(r+kr,c+kc,i); + const double w = kernel(kr,kc)*gauss(val, lval, sigma_r, ndims); + for (octave_idx_type i = 0; i < num_planes; i++) sum[i] += w*lval[i]; + k += w; + } + } + for (octave_idx_type i = 0; i < num_planes; i++) out(r,c,i) = sum[i]/k; + OCTAVE_QUIT; + } + } + + return octave_value(out); +} + +DEFUN_DLD(__bilateral__, args, , "\ +-*- texinfo -*-\n\ +@deftypefn {Loadable Function} __bilateral__(@var{im}, @var{sigma_d}, @var{sigma_r})\n\ +Performs Gaussian bilateral filtering in the image @var{im}. @var{sigma_d} is the\n\ +spread of the Gaussian used as closenes function, and @var{sigma_r} is the spread\n\ +of Gaussian used as similarity function. This function is internal and should NOT\n\ +be called directly. Instead use @code{imsmooth}.\n\ +@end deftypefn\n\ +") +{ + octave_value_list retval; + if (args.length() != 3) + { + print_usage (); + return retval; + } + + const octave_idx_type ndims = args(0).ndims(); + if (ndims != 2 && ndims != 3) + { + error("__bilateral__: only 2 and 3 dimensional is supported"); + return retval; + } + const double sigma_d = args(1).scalar_value(); + const double sigma_r = args(2).scalar_value(); + if (error_state) + { + error("__bilateral__: invalid input"); + return retval; + } + + // Take action depending on input type + if (args(0).is_real_matrix()) + { + const NDArray im = args(0).array_value(); + retval = bilateral<NDArray>(im, sigma_d, sigma_r); + } + else if (args(0).is_int8_type()) + { + const int8NDArray im = args(0).int8_array_value(); + retval = bilateral<int8NDArray>(im, sigma_d, sigma_r); + } + else if (args(0).is_int16_type()) + { + const int16NDArray im = args(0).int16_array_value(); + retval = bilateral<int16NDArray>(im, sigma_d, sigma_r); + } + else if (args(0).is_int32_type()) + { + const int32NDArray im = args(0).int32_array_value(); + retval = bilateral<int32NDArray>(im, sigma_d, sigma_r); + } + else if (args(0).is_int64_type()) + { + const int64NDArray im = args(0).int64_array_value(); + retval = bilateral<int64NDArray>(im, sigma_d, sigma_r); + } + else if (args(0).is_uint8_type()) + { + const uint8NDArray im = args(0).uint8_array_value(); + retval = bilateral<uint8NDArray>(im, sigma_d, sigma_r); + } + else if (args(0).is_uint16_type()) + { + const uint16NDArray im = args(0).uint16_array_value(); + retval = bilateral<uint16NDArray>(im, sigma_d, sigma_r); + } + else if (args(0).is_uint32_type()) + { + const uint32NDArray im = args(0).uint32_array_value(); + retval = bilateral<uint32NDArray>(im, sigma_d, sigma_r); + } + else if (args(0).is_uint64_type()) + { + const uint64NDArray im = args(0).uint64_array_value(); + retval = bilateral<uint64NDArray>(im, sigma_d, sigma_r); + } + else + { + error("__bilateral__: first input should be a real or integer array"); + return retval; + } + + return retval; +} This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <ha...@us...> - 2008-04-12 15:46:55
|
Revision: 4935 http://octave.svn.sourceforge.net/octave/?rev=4935&view=rev Author: hauberg Date: 2008-04-12 08:46:59 -0700 (Sat, 12 Apr 2008) Log Message: ----------- Added support for custom Gaussian smoothing Modified Paths: -------------- trunk/octave-forge/main/image/inst/imsmooth.m trunk/octave-forge/main/image/src/Makefile Added Paths: ----------- trunk/octave-forge/main/image/src/__custom_gaussian_smoothing__.cc Modified: trunk/octave-forge/main/image/inst/imsmooth.m =================================================================== --- trunk/octave-forge/main/image/inst/imsmooth.m 2008-04-12 13:33:56 UTC (rev 4934) +++ trunk/octave-forge/main/image/inst/imsmooth.m 2008-04-12 15:46:59 UTC (rev 4935) @@ -37,6 +37,8 @@ ## @itemx "Perona and Malik" ## @itemx "P&M" ## Smoothing using anisotropic diffusion as described by Perona and Malik. +## @item "Custom Gaussian" +## Gaussian smoothing with a spatially varying covariance matrix. ## @end table ## ## In all algorithms the computation is done in double precision floating point @@ -146,6 +148,29 @@ ## @var{J} = imsmooth(@var{I}, "p&m", 25, 0.25, @var{g}); ## @end example ## +## @strong{Custom Gaussian Smoothing} +## +## The image is smoothed using a Gaussian filter with a spatially varying covariance +## matrix. The third and fourth input arguments contain the Eigenvalues of the +## covariance matrix, while the fifth contains the rotation of the Gaussian. +## These arguments can be matrices of the same size as the input image, or scalars. +## In the last case the scalar is used in all pixels. If the rotation is not given +## it defaults to zero. +## +## The following example shows how to increase the size of an Gaussian +## filter, such that it is small near the upper right corner of the image, and +## large near the lower left corner. +## +## @example +## [@var{lambda1}, @var{lambda2}] = meshgrid (linspace (0, 25, columns (@var{I})), linspace (0, 25, rows (@var{I}))); +## @var{J} = imsmooth (@var{I}, "Custom Gaussian", @var{lambda1}, @var{lambda2}); +## @end example +## +## The implementation uses an elliptic filter, where only neighbouring pixels +## with a Mahalanobis' distance to the current pixel that is less than 3 are +## used to compute the response. The response is computed using double precision +## floating points, but the result is of the same class as the input image. +## ## @strong{References} ## ## [1] P. Perona and J. Malik, @@ -162,7 +187,7 @@ ## TODO: Implement Joachim Weickert's anisotropic diffusion (it's soo cool) -function J = imsmooth(I, name = "gaussian", varargin) +function J = imsmooth(I, name = "Gaussian", varargin) ## Check inputs if (nargin == 0) print_usage(); @@ -174,6 +199,10 @@ if ((imchannels != 1 && imchannels != 3) || tmp != 1) error("imsmooth: first input argument must be an image"); endif + if (nargin == 2 && isscalar (name)) + varargin {1} = name; + name = "Gaussian"; + endif if (!ischar(name)) error("imsmooth: second input must be a string"); endif @@ -356,6 +385,38 @@ J(:,:,i) = pm(I(:,:,i), iter, lambda, method); endfor + ##################################### + ### Custom Gaussian Smoothing ### + ##################################### + case "custom gaussian" + ## Check input + if (length (varargin) < 2) + error ("imsmooth: not enough input arguments"); + elseif (length (varargin) == 2) + varargin {3} = 0; # default theta value + endif + arg_names = {"third", "fourth", "fifth"}; + for k = 1:3 + if (isscalar (varargin {k})) + varargin {k} = repmat (varargin {k}, imrows, imcols); + elseif (ismatrix (varargin {k}) && ndims (varargin {k}) == 2) + if (rows (varargin {k}) != imrows || columns (varargin {k}) != imcols) + error (["imsmooth: %s input argument must have same number of rows " + "and columns as the input image"], arg_names {k}); + endif + else + error ("imsmooth: %s input argument must be a scalar or a matrix", arg_names {k}); + endif + if (!strcmp (class (varargin {k}), "double")) + error ("imsmooth: %s input argument must be of class 'double'", arg_names {k}); + endif + endfor + + ## Perform the smoothing + for i = imchannels:-1:1 + J(:,:,i) = __custom_gaussian_smoothing__ (I(:,:,i), varargin {:}); + endfor + ###################################### ### Mean Shift Based Smoothing ### ###################################### Modified: trunk/octave-forge/main/image/src/Makefile =================================================================== --- trunk/octave-forge/main/image/src/Makefile 2008-04-12 13:33:56 UTC (rev 4934) +++ trunk/octave-forge/main/image/src/Makefile 2008-04-12 15:46:59 UTC (rev 4935) @@ -12,8 +12,9 @@ IMAGEMAGICK=__magick_read__.oct endif -all: __cordfltn__.oct __bilateral__.oct bwlabel.oct bwfill.oct rotate_scale.oct \ - hough_line.oct graycomatrix.oct deriche.oct __bwdist.oct nonmax_supress.oct \ +all: __cordfltn__.oct __bilateral__.oct __custom_gaussian_smoothing__.oct \ + bwlabel.oct bwfill.oct rotate_scale.oct hough_line.oct \ + graycomatrix.oct deriche.oct __bwdist.oct nonmax_supress.oct \ $(JPEG) $(PNG) $(IMAGEMAGICK) jpgread.oct: jpgread.cc Added: trunk/octave-forge/main/image/src/__custom_gaussian_smoothing__.cc =================================================================== --- trunk/octave-forge/main/image/src/__custom_gaussian_smoothing__.cc (rev 0) +++ trunk/octave-forge/main/image/src/__custom_gaussian_smoothing__.cc 2008-04-12 15:46:59 UTC (rev 4935) @@ -0,0 +1,204 @@ +#include <octave/oct.h> + +template <class MT> +MT +custom_gaussian_smoothing (const MT &I, const Matrix &lambda1, const Matrix &lambda2, + const Matrix &theta) +{ + const octave_idx_type rows = I.rows (); + const octave_idx_type cols = I.columns (); + + // Allocate output + MT J (I.dims ()); + + // Iterate over every element of 'I' + for (octave_idx_type row = 0; row < rows; row++) + { + for (octave_idx_type col = 0; col < cols; col++) + { + // Extract parameters + const double v1 = lambda1 (row, col); + const double v2 = lambda2 (row, col); + const double t = theta (row, col); + + // Should we perform any filtering? + if (std::min (v1, v2) > 0) + { + // Compute inverse covariance matrix, C^-1 = [a, b; b, c] + const double iv1 = 1.0/v1; + const double iv2 = 1.0/v2; + const double ct = cos (t); + const double st = sin (t); + const double ct2 = ct*ct; + const double st2 = st*st; + const double ctst = ct*st; + const double a = ct2*iv2 + st2*iv1; + const double b = (iv2-iv1)*ctst; + const double c = st2*iv2 + ct2*iv1; + + // Compute bounding box of the filter + const double k = 3.0; // The maximally allowed Mahalanobis' distance + const double sqrtv1 = sqrt (v1); + const double sqrtv2 = sqrt (v2); + const octave_idx_type rur = abs (k*(ct*sqrtv2 - st*sqrtv1)); // 'rur' means 'row-upper-right' + const octave_idx_type cur = abs (k*(st*sqrtv2 + ct*sqrtv1)); + const octave_idx_type rlr = abs (k*(ct*sqrtv2 + st*sqrtv1)); + const octave_idx_type clr = abs (k*(st*sqrtv2 - ct*sqrtv1)); + const octave_idx_type rul = abs (k*(-ct*sqrtv2 - st*sqrtv1)); + const octave_idx_type cul = abs (k*(-st*sqrtv2 + ct*sqrtv1)); + const octave_idx_type rll = abs (k*(-ct*sqrtv2 + st*sqrtv1)); + const octave_idx_type cll = abs (k*(-st*sqrtv2 - ct*sqrtv1)); + const octave_idx_type r_delta = std::max (std::max (rur, rlr), std::max (rul, rll)); + const octave_idx_type c_delta = std::max (std::max (cur, clr), std::max (cul, cll));; + // The bounding box is now (row-r_delta):(row+r_delta)x(col-c_delta):(col+c_delta). + // We, however, represent the bounding box in a local coordinate system around (row, col). + const octave_idx_type r1 = std::max (row-r_delta, 0) - row; + const octave_idx_type r2 = std::min (row+r_delta, rows-1) - row; + const octave_idx_type c1 = std::max (col-c_delta, 0) - col; + const octave_idx_type c2 = std::min (col+c_delta, cols-1) - col; + + // Perform the actual filtering + double sum = 0; + double wsum = 0; // for normalisation + for (octave_idx_type rl = r1; rl <= r2; rl++) + { + for (octave_idx_type cl = c1; cl <= c2; cl++) + { + // Compute Mahalanobis' distance + const double dsquare = rl*(a*rl + b*cl) + cl*(b*rl + c*cl); + + // We only do the filtering in an elliptical window + if (dsquare > k*k) continue; + + // Update filter values + const double w = exp (-0.5*dsquare); + wsum += w; + sum += w*(double)I.elem (row + rl, col + cl); + } // End: cl + } // End: rl + + // Compute final result + J (row, col) = sum/wsum; + } + else // No filtering is performed + { + J.elem (row, col) = I.elem (row, col); + } + } // End: column iteration + } // End: row iteration + + // Return + return J; +} + +#define RETURN_IF_INVALID \ + if (error_state) \ + { \ + error ("__custom_gaussian_smoothing__: invalid input"); \ + return retval; \ + } + +DEFUN_DLD (__custom_gaussian_smoothing__, args, ,"\ +-*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{J} =} __custom_gaussian_smooting__ (@var{I}, @var{lambda1}, @var{lambda2}, @var{theta})\n\ +Performs Gaussian smoothing on the image @var{I}. In pixel @math{(r,c)} the \n\ +Eigenvalues of the Gaussian is @var{lambda1}@math{(r,c)} and @var{lambda2}@math{(r,c)}.\n\ +The Gaussian is rotated with the angle given in @var{theta}@math{(r,c)}.\n\ +\n\ +@strong{Warning:} this function should @i{never} be called directly! The user\n\ +interface to this function is available in @code{imsmooth}.\n\ +@seealso{imsmooth}\n\ +@end deftypefn\n\ +") +{ + // Handle Input + octave_value_list retval; + const int nargin = args.length (); + + if (nargin != 4) + { + print_usage (); + return retval; + } + + const Matrix lambda1 = args (1).matrix_value (); + const Matrix lambda2 = args (2).matrix_value (); + const Matrix theta = args (3).matrix_value (); + + RETURN_IF_INVALID; + + const octave_idx_type rows = args (0).rows(); + const octave_idx_type cols = args (0).columns(); + if (lambda1.rows () != rows || lambda1.columns () != cols + || lambda2.rows () != rows || lambda2.columns () != cols + || theta.rows () != rows || theta.columns () != cols) + { + error ("__custom_gaussian_smoothing__: size mismatch"); + return retval; + } + + // Take action depending on input type + //octave_value J; + if (args(0).is_real_matrix()) + { + const Matrix I = args(0).matrix_value(); + RETURN_IF_INVALID; + retval.append (custom_gaussian_smoothing<Matrix>(I, lambda1, lambda2, theta)); + } + else if (args(0).is_int8_type()) + { + const int8NDArray I = args(0).int8_array_value(); + RETURN_IF_INVALID; + retval.append (custom_gaussian_smoothing<int8NDArray>(I, lambda1, lambda2, theta)); + } + else if (args(0).is_int16_type()) + { + const int16NDArray I = args(0).int16_array_value(); + RETURN_IF_INVALID; + retval.append (custom_gaussian_smoothing<int16NDArray>(I, lambda1, lambda2, theta)); + } + else if (args(0).is_int32_type()) + { + const int32NDArray I = args(0).int32_array_value(); + RETURN_IF_INVALID; + retval.append (custom_gaussian_smoothing<int32NDArray>(I, lambda1, lambda2, theta)); + } + else if (args(0).is_int64_type()) + { + const int64NDArray I = args(0).int64_array_value(); + RETURN_IF_INVALID; + retval.append (custom_gaussian_smoothing<int64NDArray>(I, lambda1, lambda2, theta)); + } + else if (args(0).is_uint8_type()) + { + const uint8NDArray I = args(0).uint8_array_value(); + RETURN_IF_INVALID; + retval.append (custom_gaussian_smoothing<uint8NDArray>(I, lambda1, lambda2, theta)); + } + else if (args(0).is_uint16_type()) + { + const uint16NDArray I = args(0).uint16_array_value(); + RETURN_IF_INVALID; + retval.append (custom_gaussian_smoothing<uint16NDArray>(I, lambda1, lambda2, theta)); + } + else if (args(0).is_uint32_type()) + { + const uint32NDArray I = args(0).uint32_array_value(); + RETURN_IF_INVALID; + retval.append (custom_gaussian_smoothing<uint32NDArray>(I, lambda1, lambda2, theta)); + } + else if (args(0).is_uint64_type()) + { + const uint64NDArray I = args(0).uint64_array_value(); + RETURN_IF_INVALID; + retval.append (custom_gaussian_smoothing<uint64NDArray>(I, lambda1, lambda2, theta)); + } + else + { + error("__custom_gaussian_smoothing__: first input should be a real or integer array"); + return retval; + } + + // Return + return retval; +} This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |