# HG changeset patch # User Joonas Lipping # Date 1393435056 -7200 # Wed Feb 26 19:17:36 2014 +0200 # Node ID 28b9b1bfb6e9cf98cb066d632a890c5f1f338e56 # Parent 714073ad045f0abd8a2adf50abc608dc8ed026ad Changed fspecial('disk') to mimic Matlab behavior. (bug #41186) diff -r 714073ad045f -r 28b9b1bfb6e9 inst/fspecial.m --- a/inst/fspecial.m Tue Feb 18 23:10:42 2014 +0000 +++ b/inst/fspecial.m Wed Feb 26 19:17:36 2014 +0200 @@ -119,12 +119,65 @@ else radius = 5; endif + rsq = radius*radius; ## Create the filter - [x, y] = meshgrid (-radius:radius, -radius:radius); - r = sqrt (x.^2 + y.^2); - f = (r <= radius); - ## Normalize the filter to integral 1 - f = f / sum (f (:)); + if (radius==0) + f=1; + else + ## Area between the x axis and curve x^2 + y^2 = r^2, as x = a .. b + function auc = auc (a, b) + aroot = sqrt (rsq - a*a); + aint = a*aroot + rsq*atan (a/aroot); + broot = sqrt (rsq - b*b); + bint = b*broot + rsq*atan (b/broot); + auc = (bint-aint)/2; + endfunction + [x, y] = meshgrid(-radius:radius, -radius:radius); + rhi = (abs(x)+0.5).^2 + (abs(y)+0.5).^2; + f = (rhi <= rsq) / 1.0; + ## Place the axial caps as a special case + cap = auc (-.5,.5) - (radius-0.5); + f(1,radius+1) = cap; + f(2*radius+1,radius+1) = cap; + f(radius+1,1) = cap; + f(radius+1,2*radius+1) = cap; + ## Trace a path along the rim + x=1; + y=radius; + while (x0) + if (rsq - (x-0.5)^2 < (y+0.5)^2) + a = x-0.5; + if (rsq - (x+0.5)^2 > (y-0.5)^2) + b = x+0.5; + else + b = sqrt (rsq - (y-0.5)^2); + endif + val = auc (a,b) - (b-a)*(y-0.5); + else + a = sqrt (rsq-(y+0.5)^2); + if (rsq - (x+0.5)^2 > (y-0.5)^2) + b = x+0.5; + else + b = sqrt (rsq-(y-0.5)^2); + endif + val = auc (a,b) - (b-a)*(y-0.5) + (a - (x-0.5)); + endif + f(radius+1+x,radius+1+y) = val; + f(radius+1-x,radius+1+y) = val; + f(radius+1+x,radius+1-y) = val; + f(radius+1-x,radius+1-y) = val; + if((x-0.5)^2 + (y-1.5)^2 < rsq && (x+0.5)^2 + (y-0.5)^2 > rsq) + y=y-1; + elseif((x+0.5)^2 + (y-0.5)^2 < rsq && (x+1.5)^2 + (y+0.5)^2 > rsq) + x=x+1; + else + x=x+1; + y=y-1; + endif + endwhile + ## Normalize + f = f / pi / rsq; + endif case "gaussian" ## Get hsize