BM3D/BM3D-SAPCA/function_LPAKernelMatrixTheta.m

   1 % Return the discrete kernels for LPA estimation and their degrees matrix\r
   2 %\r
   3 % function [G, G1, index_polynomials]=function_LPAKernelMatrixTheta(h2,h1,window_type,sig_wind,TYPE,theta, m)\r
   4 %\r
   5 %\r
   6 % Outputs:\r
   7 %\r
   8 % G  kernel for function estimation\r
   9 % G1 kernels for function and derivative estimation\r
  10 %       G1(:,:,j), j=1 for function estimation, j=2 for d/dx, j=3 for d/dy,\r
  11 %                contains 0 and all higher order kernels (sorted by degree:\r
  12 %                 1 x y y^2 x^3 x^2y xy^2 y^3 etc...)\r
  13 % index_polynomials  matrix of degrees first column x powers, second\r
  14 %                    column y powers, third column total degree\r
  15 %\r
  16 %\r
  17 % Inputs:\r
  18 %\r
  19 % h2, h1  size of the kernel (size of the "asymmetrical portion")\r
  20 % m=[m(1) m(2)] the vector order of the LPA    any order combination should work \r
  21 % "theta" is an angle of the directrd window\r
  22 % "TYPE" is a type of the window support\r
  23 % "sig_wind" - vector - sigma parameters of the Gaussian wind\r
  24 % "beta"- parameter of the power in some weights for the window function\r
  25 %           (these last 3 parameters are fed into function_Window2D function)\r
  26 %\r
  27 %\r
  28 % Alessandro Foi, 6 march 2004\r
  29 \r
  30 function [G, G1, index_polynomials]=function_LPAKernelMatrixTheta(h2,h1,window_type,sig_wind,TYPE,theta, m)\r
  31 global beta\r
  32 \r
  33 %G1=0;\r
  34 m(1)=min(h1,m(1));\r
  35 m(2)=min(h2,m(2));\r
  36 \r
  37 % builds ordered matrix of the monomes powers\r
  38 number_of_polynomials=(min(m)+1)*(max(m)-min(m)+1)+(min(m)+1)*min(m)/2;   %   =size(index_polynomials,1) \r
  39 index_polynomials=zeros(number_of_polynomials,2);\r
  40 index3=1;\r
  41 for index1=1:min(m)+1\r
  42     for index2=1:max(m)+2-index1\r
  43         index_polynomials(index3,:)=[index1-1,index2-1];\r
  44         index3=index3+1;\r
  45     end\r
  46 end\r
  47 if m(1)>m(2)\r
  48     index_polynomials=fliplr(index_polynomials);\r
  49 end\r
  50 index_polynomials(:,3)=index_polynomials(:,1)+index_polynomials(:,2);    %calculates degrees of polynomials\r
  51 index_polynomials=sortrows(sortrows(index_polynomials,2),3);             %sorts polynomials by degree (x first)\r
  52 \r
  53 %=====================================================================================================================================\r
  54 \r
  55 halfH=max(h1,h2);\r
  56 H=-halfH+1:halfH-1;\r
  57 \r
  58 \r
  59 % creates window function and then rotates it\r
  60 %   win_fun=zeros(halfH-1,halfH-1);\r
  61 for x1=H\r
  62     for x2=H\r
  63         if TYPE==00|TYPE==200|TYPE==300 % SYMMETRIC WINDOW\r
  64             win_fun1(x2+halfH,x1+halfH)=function_Window2D(x1/h1/(1-1000*eps),x2/h2/(1-1000*eps),window_type,sig_wind,beta,h2/h1); % weight        \r
  65         end\r
  66         if TYPE==11|TYPE==211|TYPE==311 % NONSYMMETRIC ON X1,X2 WINDOW\r
  67             win_fun1(x2+halfH,x1+halfH)=(x1>=-0.05)*(x2>=-0.05)*function_Window2D(x1/h1/(1-1000*eps),x2/h2/(1-1000*eps),window_type,sig_wind,beta,h2/h1); % weight\r
  68         end\r
  69         if TYPE==10|TYPE==210|TYPE==310 % NONSYMMETRIC ON X1 WINDOW\r
  70             win_fun1(x2+halfH,x1+halfH)=(x1>=-0.05)*function_Window2D(x1/h1/(1-1000*eps),x2/h2/(1-1000*eps),window_type,sig_wind,beta,h2/h1); % weight\r
  71         end\r
  72         \r
  73         if TYPE==100|TYPE==110|TYPE==111 % exact sampling\r
  74             xt1=x1*cos(-theta)+x2*sin(-theta);\r
  75             xt2=x2*cos(-theta)-x1*sin(-theta);\r
  76             if TYPE==100 % SYMMETRIC WINDOW\r
  77                 win_fun1(x2+halfH,x1+halfH)=function_Window2D(xt1/h1/(1-1000*eps),xt2/h2/(1-1000*eps),window_type,sig_wind,beta,h2/h1); % weight        \r
  78             end\r
  79             if TYPE==111 % NONSYMMETRIC ON X1,X2 WINDOW\r
  80                 \r
  81                 win_fun1(x2+halfH,x1+halfH)=(xt1>=-0.05)*(xt2>=-0.05)*function_Window2D(xt1/h1/(1-1000*eps),xt2/h2/(1-1000*eps),window_type,sig_wind,beta,h2/h1); % weight\r
  82             end\r
  83             if TYPE==110 % NONSYMMETRIC ON X1 WINDOW\r
  84                 \r
  85                 win_fun1(x2+halfH,x1+halfH)=(xt1>=-0.05)*function_Window2D(xt1/h1/(1-1000*eps),xt2/h2/(1-1000*eps),window_type,sig_wind,beta,h2/h1); % weight\r
  86             end\r
  87         end\r
  88         \r
  89     end\r
  90 end\r
  91 win_fun=win_fun1;\r
  92 if (theta~=0)&(TYPE<100)\r
  93     win_fun=imrotate(win_fun1,theta*180/pi,'nearest');     % use 'nearest' or 'bilinear' for different interpolation schemes ('bicubic'...?)\r
  94 end\r
  95 if (theta~=0)&(TYPE>=200)&(TYPE<300)\r
  96     win_fun=imrotate(win_fun1,theta*180/pi,'bilinear');     % use 'nearest' or 'bilinear' for different interpolation schemes ('bicubic'...?)\r
  97 end\r
  98 if (theta~=0)&(TYPE>=300)\r
  99     win_fun=imrotate(win_fun1,theta*180/pi,'bicubic');     % use 'nearest' or 'bilinear' for different interpolation schemes ('bicubic'...?)\r
 100 end\r
 101 \r
 102 \r
 103 \r
 104 % make the weight support a square\r
 105 win_fun2=zeros(max(size(win_fun)));\r
 106 win_fun2((max(size(win_fun))-size(win_fun,1))/2+1:max(size(win_fun))-((max(size(win_fun))-size(win_fun,1))/2),(max(size(win_fun))-size(win_fun,2))/2+1:max(size(win_fun))-((max(size(win_fun))-size(win_fun,2))/2))=win_fun;\r
 107 win_fun=win_fun2;\r
 108 \r
 109 \r
 110 %=====================================================================================================================================\r
 111 \r
 112 \r
 113 %%%%  rotated coordinates\r
 114 H=-(size(win_fun,1)-1)/2:(size(win_fun,1)-1)/2;\r
 115 halfH=(size(win_fun,1)+1)/2;\r
 116 h_radious=halfH;\r
 117 Hcos=H*cos(theta); Hsin=H*sin(theta);\r
 118 \r
 119 \r
 120 %%%% Calculation of FI matrix\r
 121 FI=zeros(number_of_polynomials);\r
 122 i1=0;\r
 123 for s1=H\r
 124     i1=i1+1;\r
 125     i2=0;\r
 126     for s2=H\r
 127         i2=i2+1;\r
 128         x1=Hcos(s1+h_radious)-Hsin(s2+h_radious);\r
 129         x2=Hsin(s1+h_radious)+Hcos(s2+h_radious);\r
 130         phi=sqrt(win_fun(s2+halfH,s1+halfH))*(prod(((ones(number_of_polynomials,1)*[x1 x2]).^index_polynomials(:,1:2)),2)./prod(gamma(index_polynomials(:,1:2)+1),2).*(-ones(number_of_polynomials,1)).^index_polynomials(:,3));\r
 131         FI=FI+phi*phi';\r
 132     end % end of s2\r
 133 end % end of s1\r
 134 \r
 135 %FI_inv=((FI+1*eps*eye(size(FI)))^(-1));  % invert FI matrix\r
 136 FI_inv=pinv(FI);   % invert FI matrix (using pseudoinverse)\r
 137 G1=zeros([size(H,2) size(H,2) number_of_polynomials]);\r
 138 \r
 139 %%%% Calculation of mask\r
 140 i1=0;\r
 141 for s1=H\r
 142     i1=i1+1;\r
 143     i2=0;\r
 144     for s2=H\r
 145         i2=i2+1;\r
 146         x1=Hcos(s1+h_radious)-Hsin(s2+h_radious);\r
 147         x2=Hsin(s1+h_radious)+Hcos(s2+h_radious);\r
 148         phi=FI_inv*win_fun(s2+halfH,s1+halfH)*(prod(((ones(number_of_polynomials,1)*[x1 x2]).^index_polynomials(:,1:2)),2)./prod(gamma(index_polynomials(:,1:2)+1),2).*(-ones(number_of_polynomials,1)).^index_polynomials(:,3));\r
 149         G(i2,i1,1)=phi(1);              %   Function Est\r
 150         G1(i2,i1,:)=phi(:)';  % Function est & Der est on X Y etc...\r
 151     end % end of s1\r
 152 end % end of s2\r
 153 %keyboard