diff -r 714073ad045f NEWS --- a/NEWS Tue Feb 18 23:10:42 2014 +0000 +++ b/NEWS Wed Mar 05 23:04:58 2014 +0200 @@ -1,3 +1,7 @@ + ** The 'disk' filter of the function fspecial has been changed in a + backwards incompatible way. The elements on the rim of the disk are now + weighted by how much of them is covered by the disk. + Summary of important user-visible changes for image 2.4.0 (yyyy/mm/dd): ------------------------------------------------------------------------- diff -r 714073ad045f inst/fspecial.m --- a/inst/fspecial.m Tue Feb 18 23:10:42 2014 +0000 +++ b/inst/fspecial.m Wed Mar 05 23:04:58 2014 +0200 @@ -27,8 +27,11 @@ ## created. ## @item "disk" ## Circular averaging filter. The optional argument @var{arg1} controls the -## radius of the filter. If @var{arg1} is an integer @var{N}, a 2 @var{N} + 1 -## filter is created. By default a radius of 5 is used. +## radius of the filter. If @var{arg1} is an integer @var{R}, a 2 @var{R} + 1 +## filter is created. By default a radius of 5 is used. If the returned matrix +## corresponds to a cartesian grid, each element of the matrix is weighted by +## how much of the corresponding grid square is covered by a disk of radius +## @var{R} and centered at the middle of the element @var{R}+1,@var{R}+1. ## @item "gaussian" ## Gaussian filter. The optional argument @var{arg1} controls the size of the ## filter. If @var{arg1} is an integer @var{N}, a @var{N} by @var{N} @@ -115,16 +118,90 @@ case "disk" ## Get the radius if (nargin > 1 && isreal (arg1) && length (arg1 (:)) == 1) - radius = arg1; + r = arg1; else - radius = 5; + r = 5; endif ## 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 (r==0) + f=1; + else + ax=r+1; # index of the "x-axis" and "y-axis" + corner=floor (r/sqrt (2)+0.5)-0.5; # corner corresponding to 45 degrees + rsq=r*r; + ## First set values for points completely covered by the disk + [X, Y] = meshgrid (-r:r, -r:r); + rhi = (abs (X)+0.5).^2 + (abs (Y)+0.5).^2; + f = (rhi <= rsq) / 1.0; + xx=linspace (0.5,r-0.5,r); + ii=sqrt (rsq-xx.^2); # intersection points for sqrt (r^2 - x^2) + ## Set the values at the axis caps + tmp=sqrt (rsq-0.25); + rint=(0.5*tmp+rsq*atan (0.5/tmp))/2; # value of integral on the right + cap=2*rint-r+0.5; # at the caps, lint = rint + f(ax,ax+r)=cap; + f(ax,ax-r)=cap; + f(ax+r,ax)=cap; + f(ax-r,ax)=cap; + if (r==1) + y=ii(1); + lint=rint; + tmp=sqrt (rsq-y^2); + rint=(y*tmp+rsq*atan (y/tmp))/2; + val=rint-lint-0.5*(y-0.5); + f(ax-r,ax-r)=val; + f(ax+r,ax-r)=val; + f(ax-r,ax+r)=val; + f(ax+r,ax+r)=val; + else + ## Set the values elsewhere on the rim + idx=1; # index in the vector ii + x=0.5; # bottom left corner of the current square + y=r-0.5; + ybreak=false; # did we change our y last time + rx=0.5; # x on the right of the integrable region + do + i=x+0.5; + j=y+0.5; + lint=rint; + lx=rx; + if (ybreak) + ybreak=false; + val=lx-x; + idx++; + x++; + rx=x; + val-=y*(x-lx); + elseif (ii(idx+1) < y) + ybreak=true; + y--; + rx=ii(y+1.5); + val=(y+1)*(x-rx); + else + val=-y; + idx++; + x++; + rx=x; + if (floor (ii(idx)-0.5) == y) + y++; + endif + endif + tmp=sqrt (rsq-rx*rx); + rint=(rx*tmp+rsq*atan (rx/tmp))/2; + val+=rint-lint; + f(ax+i, ax+j)=val; + f(ax+i, ax-j)=val; + f(ax-i, ax+j)=val; + f(ax-i, ax-j)=val; + f(ax+j, ax+i)=val; + f(ax+j, ax-i)=val; + f(ax-j, ax+i)=val; + f(ax-j, ax-i)=val; + until (ycorner) + endif + # Normalize + f /= pi*rsq; + endif case "gaussian" ## Get hsize @@ -250,3 +327,35 @@ error ("fspecial: filter type '%s' is not supported", type); endswitch endfunction + + +## Test that the disk filter's error does not grow unreasonably large +%!test +%! eps=1e-15; +%! for i=1:9 +%! n=2^i; +%! assert(sum (fspecial ('disk', n)(:)), 1, eps*n*n); +%! endfor + +## Test that all squares completely under the disk or completely out of it are +## being assigned the correct values. +%!test +%! eps=1e-15; +%! for r=[3,5,9,17] +%! f = fspecial ('disk', r); +%! [X, Y] = meshgrid (-r:r, -r:r); +%! rhi = (abs (X)+0.5).^2 + (abs (Y)+0.5).^2; +%! rlo = (abs (X)-0.5).^2 + (abs (Y)-0.5).^2; +%! fhi = (rhi <= (r^2)); +%! flo = (rlo >= (r^2)); +%! for i=1:(2*r+1) +%! for j=1:(2*r+1) +%! if (fhi(i,j)) +%! assert(f(i,j), 1/(pi*r^2), eps); +%! endif +%! if (flo(i,j)) +%! assert(f(i,j), 0); +%! endif +%! endfor +%! endfor +%! endfor