fspecial滤波算子

Halcom

1. function h = fspecial(type,p2)
   
2. 
   
3. if nargin < 2
   
4.     switch type
   
5.         case 'average'
   
6.         p2 = [3 3];  % siz
   
7.     case 'disk'
   
8.         p2 = 5;      % rad
   
9.     end
   
10. end
    
11. 
    
12. switch type
    
13.   case 'average' % Smoothing filter
    
14.      siz = p2;
    
15.      h   = ones(siz)/prod(siz);
    
16. 
    
17.   case 'disk' % Disk filter
    
18.      rad   = p2;
    
19.      crad  = ceil(rad-0.5);
    
20.      [x,y] = meshgrid(-crad:crad,-crad:crad);
    
21.      maxxy = max(abs(x),abs(y));
    
22.      minxy = min(abs(x),abs(y));
    
23.      m1 = (rad^2 <  (maxxy+0.5).^2 + (minxy-0.5).^2).*(minxy-0.5) + ...
    
24.           (rad^2 >= (maxxy+0.5).^2 + (minxy-0.5).^2).* ...
    
25.                 sqrt(rad^2 - (maxxy + 0.5).^2);
    
26.      m2 = (rad^2 >  (maxxy-0.5).^2 + (minxy+0.5).^2).*(minxy+0.5) + ...
    
27.           (rad^2 <= (maxxy-0.5).^2 + (minxy+0.5).^2).* ...
    
28.            sqrt(rad^2 - (maxxy - 0.5).^2);
    
29.      sgrid = (rad^2*(0.5*(asin(m2/rad) - asin(m1/rad)) + ...
    
30.              0.25*(sin(2*asin(m2/rad)) - sin(2*asin(m1/rad)))) - ...
    
31.              (maxxy-0.5).*(m2-m1) + (m1-minxy+0.5)) ...
    
32.                   .*((((rad^2 < (maxxy+0.5).^2 + (minxy+0.5).^2) & ...
    
33.              (rad^2 > (maxxy-0.5).^2 + (minxy-0.5).^2)) | ...
    
34.                   ((minxy==0)&(maxxy-0.5 < rad)&(maxxy+0.5>=rad))));
    
35.      sgrid = sgrid + ((maxxy+0.5).^2 + (minxy+0.5).^2 < rad^2);
    
36.      sgrid(crad+1,crad+1) = min(pi*rad^2,pi/2);
    
37.      if ((crad>0) && (rad > crad-0.5) && (rad^2 < (crad-0.5)^2+0.25))
    
38.         m1  = sqrt(rad^2 - (crad - 0.5).^2);
    
39.              m1n = m1/rad;
    
40.         sg0 = 2*(rad^2*(0.5*asin(m1n) + 0.25*sin(2*asin(m1n)))-m1*(crad-0.5));
    
41.         sgrid(2*crad+1,crad+1) = sg0;
    
42.         sgrid(crad+1,2*crad+1) = sg0;
    
43.         sgrid(crad+1,1)        = sg0;
    
44.         sgrid(1,crad+1)        = sg0;
    
45.         sgrid(2*crad,crad+1)   = sgrid(2*crad,crad+1) - sg0;
    
46.         sgrid(crad+1,2*crad)   = sgrid(crad+1,2*crad) - sg0;
    
47.         sgrid(crad+1,2)        = sgrid(crad+1,2)      - sg0;
    
48.         sgrid(2,crad+1)        = sgrid(2,crad+1)      - sg0;
    
49.      end
    
50.      sgrid(crad+1,crad+1) = min(sgrid(crad+1,crad+1),1);
    
51.      h = sgrid/sum(sgrid(:));
    
52. 
    
53. end
    
54.   
    
55. end
    

复制代码