diff -r 714073ad045f NEWS --- a/NEWS Tue Feb 18 23:10:42 2014 +0000 +++ b/NEWS Sat Mar 01 23:29:31 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 Sat Mar 01 23:29:31 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,93 @@ 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) + %printf('post-ybreak\n'); + ybreak=false; + val=lx-x; + idx++; + x++; + rx=x; + val-=y*(x-lx); + elseif (ii(idx+1) < y) + %printf('ybreak\n'); + ybreak=true; + y--; + rx=ii(y+1.5); + val=(y+1)*(x-rx); + else + %printf('xbreak\n'); + 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 +330,18 @@ 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; +%! err(i)=abs(1-sum(fspecial('disk',n)(:))); +%! # Error is due to double precision; it is proportional to the number of +%! # elements in the matrix. +%! maxerr(i)=eps*n^2; +%! endfor +%! if (!all(err < maxerr)) +%! error("Unexpected large error in fspecial disk filter"); +%! endif