#-------------------------------------------------------------------------------------------------------------------------
# author :  zijia li  
# date   :  Jan,  2010.
#--------------------------------------------------------------------------------------------------------------------------



###########################################################################################################################
## The code is used to compute the structure matrix about the given  M1 and M2 from one RGB version by fast discrete 
## Fouriertransform where the algorithm= mintime mesn it is FFT.
## Input:
##   Pho1 ---blurred RGB version
##     e----the tolerance  
## Output:
##     Phoo----deblurring version of one RGB version
########################################################################################################################## 

  
##########################################################################################################################
## Firstly,we must get the code of approximate 1-D GCD which is useful in the 2-D GCD algorithm.Secondly,we should get the code 
## of compute matrix A and B which are very useful in computing the coefficients of 2-D GCD. Thirdly,we should get the 
## code of compute matrix TT which is used to compute the coefficients of 2-D GCD directly using svd decomposition.
## At last we can get the code of completely.
##########################################################################################################################




########################################################################################################################### 
## In Ap1D_gcd
##     in this algorithm, we get the polynomial approximate gcd of f and g utilizing SVD for computing cofactor and Dividen of 
##           DFT : 
##
## Input:
##     f,p----polynomials
##      e ----tolerance
##      deg --- the degree of gcd
##       
## Output 
##     G ----polynomials we need 
##     CC---- cofactor of f
##     DD---- cofactor of g
##      d---- the degree of gcd
###########################################################################################################################


Ap1D_gcd:=proc( f :: polynom, g :: polynom, x :: name, e :: numeric, deg :: numeric)
  local 
    m_f,n_g                                   # the degree of the polynomials
    ,SM,EE                                    # Sylvester matrix of two polynomials
    ,X,K,Y,XX,YY,H,HH,Z1,Z2,Vf,Vg             # the vectors
    ,Cv,Dv                                    # the vectors
    ,G,G1,G2,CC,DD                            # the return result and values
    ,dim                                      # the dimension of the matrix
    ,dim_zero                                 # the number of zeros
    ,_S                                       # the singular values of the matrix
    ,i,j,k,l                                  # the parameters in for loop statement 
    ,dim_row                                  # the  value in computing
    ,dim_col                                  # the  column dimension of the matrix
    ,d                                        # the degree of gcd of f and g
    ,aa,bb,cc                                 # the elements of the matrix
    ,UU,SS,VV,_U,_V                           # singular value decomposition 
    ,ff,gg                                    # the polynomials
    ,st
    ;
    ff:=expand(f/lcoeff(f,x));
    gg:=expand(g/lcoeff(g,x));
    m_f:=degree(ff);
    n_g:=degree(gg);    
    Vf:=CoefficientVector(ff,x);
    Vg:=CoefficientVector(gg,x);
    SM:=Matrix(m_f+n_g,m_f+n_g);
    for i from 1 to n_g do
      for j from 1 to m_f+1 do
        SM[i,i-1+j]:=Vf[m_f+1-j+1];
      end do;
    end do;
    for i from 1 to m_f do
      for j from 1 to n_g+1 do
        SM[n_g+i,m_f-i+j]:=Vg[n_g+1-j+1];
      end do;
    end do;
    d:=0;
    aa:=0;
    userinfo(3,S_stls,print(`construct the Sylvester matrix SM but cannot use the code SylvesterMatrix`),print('SM'=SM));
    if deg=0 then
      UU,SS,VV:=LinearAlgebra:-SingularValues(SM,output=['U','S','Vt']);
      _S:=SS;
      userinfo(3,S_stls,print(`so the singular values list is: `),print('_S'=SS));
      while _S[m_f+n_g-d]<=10^(-e) and d<min(m_f,n_g) do
        d:=d+1;
      end do;
      userinfo(3,S_stls,print(`so r (the degree of the gcd of f and g) is `),print('r'=d));
      if d=0 then
        return 1,ff/norm(expand(ff),2),gg/norm(expand(gg),2),d;
      end if;
    else
      d:=deg;
    end if;
    if d=n_g then
      return gg/norm(expand(gg),2),Deconvolution(ff,gg,x,e),1,d;
    end if;
    if d=m_f then 
      return ff/norm(expand(ff),2),1,Deconvolution(gg,ff,x,e),d;
    end if;
    XX:=Vector(m_f-d+1,i->x^(i-1));
    YY:=Vector(n_g-d+1,i->x^(i-1));
    if d=1 then 
      UU:=LinearAlgebra:-SingularValues(SM,output=['U']);
      _U:=UU;
      HH:=LinearAlgebra:-HermitianTranspose(Column(_U,m_f+n_g));
    else
      dim_row:=m_f+n_g-2*d+2;
      dim_col:=m_f+n_g-d+1;
      EE:=Matrix(dim_row,dim_col);
      for i from 1 to dim_row do
        for j from 1 to dim_col do
	  EE[i,j]:=SM[d-1+i,d-1+j];
	  end do;
      end do;
      userinfo(3,S_stls,print(`so S (the submatrix of the sylvestermatrix of f and g) is:`),print('S'=EE));
      UU:=LinearAlgebra:-SingularValues(EE,output=['U']);
      _U:=UU;
      HH:=LinearAlgebra:-HermitianTranspose(Column(_U,dim_row));
    userinfo(3,S_stls,print(`so norm of xS=0 is:`),print('Norm'=Norm(HH.EE,2)));
    end if;

    userinfo(3,S_stls,print(`construct the singular vector for the coefficients of the C(x) and D(x) which are the co-prime parts of f and g respectly`),print('H'=HH));
    DD:=0;
    for i from 1 to n_g-d+1 do
      DD:=DD+HH[i]*x^(n_g-d+1-i);
    end do;
    CC:=0;
    for i from 1 to m_f-d+1 do
      CC:=CC-HH[n_g-d+1+i]*x^(i-1);
    end do;
    CC:=expand(CC/lcoeff(CC));
    DD:=expand(DD/lcoeff(DD));
    st:=time():
    G:=Dividen_1D(f,CC,x):
    G:=G/lcoeff(G):
    st:=time()-st:
    userinfo(1,S_stls,print(`the time for gcd from cofactor`),print('t'=st));
    userinfo(3,S_stls,print(`norm :`),print('Norm'=expand(G)));
    return G/norm(expand(G),2),CC,DD,d;
end proc;
########################################################################################################################### 
## In Dividen,
##     in this algorithm, we get the polynomial CC is f/p: 
##
## Input:
##     f,p----polynomials
##       
## Output 
##     CC ----polynomials we need 
###########################################################################################################################
Dividen_1D:=proc(f,p,x)
   local
     m
     ,r
     ,i     
     ,X
     ,Vectp,Vectf_x,Vectp_x
     ,V1,V2,V3,V4
     ,CC
     ,st
     ;
     m:=degree(f,x)+1:
     r:=degree(p,x)+1:
     st:=time():
     Vectf_x:=CoefficientVector(f,x);
     Vectp:=CoefficientVector(p,x);
     Vectp_x:=Vector(m):
     for i from 1 to r do
       Vectp_x[i]:=Vectp[i];
     end do;
     X :=Vector(m, i->x^(i-1)); 
     st:=time()-st:
     userinfo(3,S_stls,print(`we get Polynomials' coefficient vector using computing time is `),print('st1'=st));
     st:=time():
     V1:=FourierTransform(Vectf_x, inplace = true, algorithm= mintime);
     V2:=FourierTransform(Vectp_x, inplace = true, algorithm= mintime);
     V3:=Vector(m);
     st:=time()-st:
     userinfo(3,S_stls,print(`we get FourierTransform using computing time is `),print('st2'=st));
     st:=time():
     for i from 1 to m do
       V3[i]:=V1[i]/V2[i];
     end do;
     V4:=DiscreteTransforms:-InverseFourierTransform(V3, inplace = true, algorithm= mintime);
     for i from m-r+2 to m do
         V4[i]:=0;
     end do;
     st:=time()-st:
     userinfo(3,S_stls,print(`we get dividen using computing time is `),print('st3'=st));     
     CC:=LinearAlgebra:-Transpose(X).V4;
     return CC;
 end proc:
########################################################################################################################### 
## In Get_gcddegree,
##     in this algorithm, we may compute the degree of the 1D degree: 
##                       
## Input:
##     f,g ----the 1-D polynomials
##       
## Output 
##     deg-----the degre of the gcd of f and g
##########################################################################################################################


Get_gcddegree:=proc(f,g,x,e)
   # e is the tolerance times
   local 
     m,n              # the degree of f and g
     ,MM              # the sylvester matrix
     ,_S               # the singularvalues of the sylvestermatrix
     ,i
     ,deg             # the dgree we need
     ;
     m:=degree(f,x);
     n:=degree(g,x);
     MM:=LinearAlgebra:-SylvesterMatrix(f,g,x);
     _S:=LinearAlgebra:-SingularValues(MM,output='list');
     deg:=0;
     userinfo(3,S_stls,print(`the interpolation and compute inverse matrix time is `),print('st1'=st),print('degx'=degx),print('degy'=MM)); 
     while _S[m+n-deg]<=10^(-e) and deg<min(m,n) do
       deg:=deg+1;
     end do:
     return deg;
   end proc:
  

########################################################################################################################### 
## In Get_AB,
##     in this algorithm, we may compute two cases(like the example in the paper): 
##     case 1: first take x by the exp(-2*Pi*k*I/m),then get the polynomials about y,using the co1D_GCD compute the 1D Gcd 
##             of F1 and F2 ,then take y by the exp(-2*Pi*l*I/m),then we get A[k,l],then in looping we get A.
##     case 2: first take y by the exp(-2*Pi*l*I/m),then get the polynomials about x,using the co1D_GCD compute the 1D Gcd 
##             of F1 and F2 ,then take x by the exp(-2*Pi*k*I/m),then we get B[k,l],then in looping we get B.
##                  
## Input:
##     M_f,M_g ----the coeff matrices of 2-D polynomials
##       
## Output 
##     A,B----matrixs which satisfy the conditions of 2D gcd
###########################################################################################################################

Get_AB:=proc( M_f,M_g, e :: numeric)
   local 
     X,Y                                            # the vectors of x
     ,A,B                                           # the matrixs we need
     ,i,j                                           # the parameters in for loop statement 
     ,m,n                                           # the rows and columns of coeff matrices M_f and M_g
     ,Z,ZZ,Z3,ZZ3                                   # the DFT matriices 
     ,F1,F2,FF1,FF2                                 # the polynomials
     ,F,FF,aaa, bbb, it                             # the return of 1D gcd
     ,degx,degy,deg                                 # the degree of 1D-GCD
     ,st,st2                                        # the time                   
     : 
     st:= time();
     m:=LinearAlgebra:-RowDimension(M_f):
     n:=LinearAlgebra:-ColumnDimension(M_f):
     Z3  :=Generate_matrixZZZZ(2,m);
     ZZ3 :=Generate_matrixZZZZ(2,n);
     X :=Vector(m, i->x^(i-1)); 
     Y :=Vector(n, i->x^(i-1));
     st:= time();
     F1 := Row(Z3, 1).M_f.Y;
     F2 := Row(Z3, 1).M_g.Y;
     F1:=F1/norm(expand(F1),2):F2:=F2/norm(expand(F2),2):
     degy:=Get_gcddegree(F1,F2,x,e);
     FF1 := LinearAlgebra:-Transpose(M_f.Column(LinearAlgebra:-Transpose(ZZ3),1)).X;
     FF2 := LinearAlgebra:-Transpose(M_g.Column(LinearAlgebra:-Transpose(ZZ3),1)).X;
     FF1:=FF1/norm(expand(FF1),2):FF2:=FF2/norm(expand(FF2),2):
     degx:=Get_gcddegree(FF1,FF2,x,e);
     st:= time()-st;
     userinfo(1,S_stls,print(`the degree of GCD is`),print('degx'=degx),print('degy'=degy),print('st'=st));
     st:= time();
     Z  := Generate_matrixZZZZ(1+degx,m);
     ZZ := Generate_matrixZZZZ(1+degy,n);
     A:=Matrix(1+degx,1+degy);
     B:=Matrix(1+degx,1+degy);
     st:= time()-st;
     userinfo(3,S_stls,print(`the interpolation vector time is `),print('st1'=Z,ZZ));
     st:= time();
     deg:=degy;
     for i from 1 to 1+degx do
        st2:= time();
       F1 := Row(Z, i).M_f.Y;
       F2 := Row(Z, i).M_g.Y;
       st2:= time()-st2;
       userinfo(1,S_stls,print(`evaluation time 1D f and g  is `),print('st1'=st2),print('deg'=degy));
       st2:=time():
       F,bbb,aaa,degy:= Ap1D_gcd(F2, F1, x, e,deg);
         if F=1 then
           for j from 1 to 1+degy do
             A[i,j] := 1;
           end do;
	 else
           for j from 1 to 1+degy do
             A[i,j] := eval(F,x=evalf(cos(2*Pi*(j-1)/(1+degy))-I*sin(2*Pi*(j-1)/(1+degy))));
           end do;
         end if;
        st2:= time()-st2;
       userinfo(1,S_stls,print(`the GCD's computing time 1D f and g  is `),print('st1'=st2),print('deg'=degy));
     end do; 
     userinfo(3,S_stls,print(`generate the matrix A :`), print('A'=A));
      st:= time()-st;
     userinfo(1,S_stls,print(`the first time GCD's computing time is `),print('st1'=st));
     st:= time();
     deg:=degx;
     for i from 1 to 1+degy do
       st2:= time();     
       FF1 := LinearAlgebra:-Transpose(X).(M_f.LinearAlgebra:-Transpose(Row(ZZ,i)));
       FF2 := LinearAlgebra:-Transpose(X).(M_g.LinearAlgebra:-Transpose(Row(ZZ,i)));
       userinfo(3,S_stls,print(`generate the 1D polynomials :`), print('F1'=FF1),print('F2'=FF2));
       st2:= time()-st2;
       userinfo(1,S_stls,print(`evaluation time 1D f and g  is `),print('st1'=st2),print('deg'=degy));
       st2:=time():       
       FF,bbb,aaa,degx := Ap1D_gcd(FF2, FF1, x, e,deg);
         if FF=1 then
           for j from 1 to 1+degx do
             B[j,i] := 1;
           end do;
         else
	   for j from 1 to 1+degx do
             B[j,i] := eval(FF,x=evalf(cos(2*Pi*(j-1)/(1+degx))-I*sin(2*Pi*(j-1)/(1+degx))));
           end do;
         end if;
        st2:= time()-st2;
       userinfo(1,S_stls,print(`the GCD's computing time 1D f and g  is `),print('st1'=st2),print('deg'=degx));
     end do; 
     st:= time()-st;
     userinfo(1,S_stls,print(`the second time GCD's computing time is `),print('st1'=st));
     userinfo(3,S_stls,print(`generate the matrix B :`), print('B'=B));
    
     return A,B;
end proc:



########################################################################################################################### 
## In Get_TT,
##     in this algorithm, we may compute the TT matrix : 
##
## Input:
##     A,B----Matrixs
##       
## Output 
##     TT ----matrix which satisfy the conditions
###########################################################################################################################
Get_TT:=proc( A , B  )
   local 
     TT                                 # the matrix
     ,i,j                               # the parameters in for loop statement 
     ,m,n                               # the rows and columns of A&B
     ,W:                                # the vector
     m:=LinearAlgebra:-RowDimension(A); 
     n:=LinearAlgebra:-ColumnDimension(A);
     TT:=Matrix(m+n,m+n);
     for i from 1 to m do 
       for j from 1 to n do 
         TT[i, i] := TT[i, i]+(abs(A[i, j]))^2; 
         TT[m+j, m+j] := TT[m+j, m+j]+(abs(B[i, j]))^2; 
         TT[i, m+j] := -conjugate(A[i, j])*B[i, j]; 
         TT[m+j, i] := -conjugate(B[i, j])*A[i, j]:                               
       end do;
     end do; 
     return TT;
end proc:

########################################################################################################################### 
## In Get_2Dgcd,
##     in this algorithm, get the coeff matrix of 2D gcd for giving two coeff matrices : 
##
## Input:
##     M1,M2---coeff matrices of 2D polynomials
##     e   ----- the tolerence
##     
##       
## Output 
##     CC ----the coeff matrix of 2D gcd we need
###########################################################################################################################

Get_2Dgcd:=proc( M1,M2, e)
   local 
     AA,BB                               # the matrixes as the defined of A & B in the paper
     ,T_T                                # the matrixes as the defined of TT in the paper
     ,Z,ZZ                               # the matrixes as the defined of TT in the paper 
     ,V1                                 # the matrix V of the Singular decomposition of TT as the defined in the paper
     ,P1,P2                              # the matrixes of P and InverseFourierTransform its as defined in the paper 
     ,m,n                                # the row and column of Matrix A&B defined in the paper
     ,i,j                                # the parameters in computing for loop statement
     ,deg,xdeg,ydeg                      # the degree of 1D-GCD
     ,X,Y                                # the vector of variable x and y 
     ,W                                  # the singular vector we need for the equation TX=0   
     ,st                                 # the time beginning
     ,CC
      :
     st:= time():
     AA,BB:=Get_AB(M1,M2,e);
     m:=LinearAlgebra:-RowDimension(AA); 
     n:=LinearAlgebra:-ColumnDimension(AA);
     P1:=Matrix(m,n);
     # get matrices A,B which are used to interpolation the CC
     st:= time():
     T_T:=Get_TT(AA,BB);
     st:= time()-st;
     userinfo(3,S_stls,print(`we get TT using computing time is `),print('st1'=st));
     st:= time():
     V1:=LinearAlgebra:-SingularValues(T_T, output = ['Vt']); 
     W:=LinearAlgebra:-HermitianTranspose(Row(V1, m+n));
     st:= time()-st;
     userinfo(3,S_stls,print(`we get TT's singularvector using computing time is `),print('st1'=st));
     st:= time():
     for i from 1 to m do
       for j from 1 to n do
         P1[i, j] := (AA[i,j]*W[i]+BB[i,j]*W[m+j])/2;
      end do; 
     end do; 
     # solve the least square
     P2:=DiscreteTransforms:-InverseFourierTransform(P1, inplace = true, algorithm= mintime);
     userinfo(3,S_stls,print(`the matrix of evaluation we need as defined in the paper`),print('P2'=P2));
     return P2;
  end proc:
########################################################################################################################### 
## In SD_Dividen,
##     in this algorithm, we get the coeff matrix of polynomial CC is f/p: 
##
## Input:
##     M_f,M_p----the coeff matrices of polynomials
##       
## Output 
##     CC ----the coeff matrix polynomial we need 
###########################################################################################################################
SD_Dividen:=proc(M_f,M_p)
   local
     m,n,mm
     ,r,s
     ,i,j
     ,X,Y
     ,Vectf_y,Vectp_y,Vectf_x,Vectp_x
     ,Z1,ZZ1,Z2,ZZ2
     ,M1,M2,M3,M4,M5
     ,CC
     ,evlu
     ,st
     ;
     st:=time():
     m:=LinearAlgebra:-RowDimension(M_f):
     n:=LinearAlgebra:-ColumnDimension(M_f):
     r:=LinearAlgebra:-RowDimension(M_p):
     s:=LinearAlgebra:-ColumnDimension(M_p):     
     X :=Vector(m, i->x^(i-1)); 
     Y :=Vector(n, i->y^(i-1));
     st:=time()-st:
    userinfo(1,S_stls,print(`we get Polynomials' coefficient matrix using computing time is `),print('st1'=st));
     st:=time():
     M1:=DiscreteTransforms:-FourierTransform(M_f, algorithm= mintime);
     st:=time()-st:
     userinfo(1,S_stls,print(`The first FFT time is `),print('st1'=st));
     st:=time():
     M5:=Matrix(m,n):
     for i from 1 to r do
       for j from 1 to s do
         M5[i,j]:=M_p[i,j];
       end do;
     end do;
     M2:=DiscreteTransforms:-FourierTransform(M5, inplace = true, algorithm= mintime);
     st:=time()-st:
     userinfo(1,S_stls,print(`The second FFT time is `),print('st1'=st));
     M3:=Matrix(m,n);
     st:=time():
     # f, p 's FourierTransform     
     for i from 1 to m do
       for j from 1 to n do
         M3[i,j]:=M1[i,j]/M2[i,j];
       end do;
     end do;
     st:=time()-st:
     userinfo(1,S_stls,print(`The divide time is `),print('st1'=st));
     st:=time():
     M4:=DiscreteTransforms:-InverseFourierTransform(M3, inplace = true, algorithm= mintime);
     st:=time()-st:
     userinfo(1,S_stls,print(`The inverse FFT time is `),print('st1'=st));
     st:=time():
     M5:=LinearAlgebra:-SubMatrix(M4,1..m-r+1,1..n-s+1);
     evlu:=abs(M5[m-r+1,n-s+1]);
     for i from 1 to m-r+1 do
       for j from 1 to n-s+1 do
         M5[i,j]:=abs(M5[i,j])/evlu;
       end do;
     end do;
     # f/p and its InverseFourierTransform
     userinfo(3,S_stls,print(`we get Polynomials using computing time is `),print('st1'=1));
     # the coeff matrix we need
     return M5;
 end proc:
########################################################################################################################### 
## In RGB_SSD_ImageDeconvolution,
##   we will deconvolution the origin version from the input two distorted version
##
##   Input  Pho1---- the blurred image versions(RGB)
##              e    ---- the accuracy
##
##   Output  Pho    ---- the version we deconvolution
########################################################################################################################### 
RGB_SSD_ImageDeconvolution:=proc(Pho1,e)
   local
     m,n                                        # the Hight and width of  version
     ,img_r1,img_g1,img_b1                      # the input Pho1's R G B
     ,f,g,h                                     # the 2D polynomials
     ,p,u1,u2,u3                                # the polynomial is 2D gcd of f and g and cofactor
     ,img_r,img_g,img_b                         # the output's R G B
     ,evlu,evlu1,evlu2,evlu3                    # the value for computing
     ,M1,M2,M3
     ,aevlu,aevlu1,aevlu2,aevlu3                # the value for computing
     ,Pho                                       # the output we need
     ,st
     ;
     
     st:=time():
     m:=ImageTools:-Height(Pho1);
     n:=ImageTools:-Width(Pho1);
     
     # the dimension of these two picture
     img_r1:= ImageTools:-Create(m, n, (i,j) -> Pho1[i,j,1]):
     img_g1:= ImageTools:-Create(m, n, (i,j) -> Pho1[i,j,2]):
     img_b1:= ImageTools:-Create(m, n, (i,j) -> Pho1[i,j,3]):
      # and there R G B to grey photo
     M1:=SD_Generate_IMmatrix(img_r1):
     M2:=SD_Generate_IMmatrix(img_g1):
     M3:=SD_Generate_IMmatrix(img_b1):
     evlu1:= RandomTools[Generate](integer(range=1..10))/10;
     evlu2:= RandomTools[Generate](integer(range=1..10))/10;
     evlu3:= RandomTools[Generate](integer(range=1..10))/10;
     aevlu1:= RandomTools[Generate](integer(range=1..10))/10;
     aevlu2:= RandomTools[Generate](integer(range=1..10))/10;
     aevlu3:= RandomTools[Generate](integer(range=1..10))/10;
     st:=time()-st:
     userinfo(1,S_stls,print(`Before we compute gcd  using computing time is `),print('st1'=st),print('evlu1'=evlu1),print('evlu2'=evlu2),print('evlu3'=evlu3),print('evlu1'=aevlu1),print('evlu2'=aevlu2),print('evlu3'=aevlu3));
     p:=Get_2Dgcd(evlu1*M1+evlu2*M2+evlu3*M3,aevlu1*M1+aevlu2*M2+aevlu3*M3,e);
     # get one of gcd
     st:=time():
     u1:=SD_Dividen(M1,p);
     u2:=SD_Dividen(M2,p);
     u3:=SD_Dividen(M3,p);
     evlu1:=Pho1[m,n,1]:
     evlu2:=Pho1[m,n,2]:
     evlu3:=Pho1[m,n,3]:
     img_r:=SD_Generate_matrixIM(u1,(6)*evlu1):
     img_g:=SD_Generate_matrixIM(u2,(6)*evlu2):
     img_b:=SD_Generate_matrixIM(u3,(6)*evlu3):
     Pho:= ImageTools:-CombineLayers(img_r, img_g, img_b):
     st:=time()-st:
     userinfo(1,S_stls,print(`AFTER we compute gcd  using computing time is `),print('st1'=st));
     return Pho;
  end proc:
########################################################################################################################### 
## In SD_Blurred,
##   we will blur the image Pho into one version from P1
##
##   Input  Pho --- the image need to be blurred
##          P1,Noise1--- the blur polynomial and noise polynomial
##
##   Output  Pho1 ----the RGB versions we need
########################################################################################################################### 
RGB_Blurred:=proc(Pho,P1,Nois1)
   local
     m,n                                       # the dimension Pho and blur filters
     ,img_r,img_g,img_b                         # the trip-pho pictures
     ,Pr,Pg,Pb                                  # the polynomials of these trip-pho
     ,Nr1,Nr2,Nr3,Ng1,Ng2,Ng3,Nb1,Nb2,Nb3       # the two versions of R G B 
     ,Pho1,Pho2,Pho3                            # the return
     ,i,j
     ;
     m:=ImageTools:-Height(Pho);
     n:=ImageTools:-Width(Pho);
     img_r := ImageTools:-Create(m, n, (i,j) -> Pho[i,j,1]):
     img_g := ImageTools:-Create(m, n, (i,j) -> Pho[i,j,2]):
     img_b := ImageTools:-Create(m, n, (i,j) -> Pho[i,j,3]):
     # the RGB versions of Pho
     Pr:=Generate_IMPolynomial(img_r,x,y):
     Nr1:=Generate_IMPolynomialmatrix(Pr*P1+Nois1,1,x,y):
     # the blurred R versions
     Pg:=Generate_IMPolynomial(img_g,x,y):
     Ng1:=Generate_IMPolynomialmatrix(Pg*P1+Nois1,1,x,y):
     # the blurred G versions
     Pb:=Generate_IMPolynomial(img_b,x,y):
     Nb1:=Generate_IMPolynomialmatrix(Pb*P1+Nois1,1,x,y):
     # the blurred B versions
     Pho1:= ImageTools:-CombineLayers(Nr1, Ng1, Nb1):
     # the combine RGB versions
   return Pho1;
   end proc:
########################################################################################################################## 
## In SNR,
##   we will compute the SNR of Pho and Noise 
##
##   Input  Pho --- the image 
##          Noise--- the noise polynomial
##
##   Output  SNR ----the values we need
########################################################################################################################### 
RGB_SNR:=proc(Pho,P1,Nois1)
   local
     m,n,m1,n1,m2,n2                # the dimension Pho and blur filters
     ,img_r,img_g,img_b             # the trip-pho pictures
     ,Pr,Pg,Pb                      # the polynomials of these trip-pho
     ,Nr1,Nr2,Nr11,Nr22       # the two versions of R G B 
     ,Pho1,Pho2                     # the return
     ,i,j
     ,evlu1,evlu2
     ,snr1r,snr2r,snr1g,snr2g,snr1b,snr2b
     ;
     m:=ImageTools:-Height(Pho);
     n:=ImageTools:-Width(Pho);
     img_r := ImageTools:-Create(m, n, (i,j) -> Pho[i,j,1]):
     img_g := ImageTools:-Create(m, n, (i,j) -> Pho[i,j,2]):
     img_b := ImageTools:-Create(m, n, (i,j) -> Pho[i,j,3]):
     # the RGB versions of Pho
     snr1r:=SD_SNR(img_r,P1,Nois1);
#     snr2r:=SD_SNR(img_r,P2,Nois2);
#     snr1g:=SD_SNR(img_g,P1,Nois1);
#     snr2g:=SD_SNR(img_g,P2,Nois2);
#     snr1b:=SD_SNR(img_b,P1,Nois1);
#     snr2b:=SD_SNR(img_b,P2,Nois2);
         
     # the combine RGB versions
   return snr1r;#snr2r,snr1g,snr2g,snr1b,snr2b;
   end proc:
########################################################################################################################### 
## In SD_SNR,
##   we will compute the SNR of Pho and Noise 
##
##   Input  Pho --- the image 
##          Noise--- the noise polynomial
##
##   Output  SNR ----the ration values we need
########################################################################################################################### 
SD_SNR:=proc(Pho,P1,Nois1)
   local
     m,n                            # the dimension Pho and blur filters
     ,Pr                            # the polynomials of these trip-pho
     ,N1,N11                        # the two versions of R G B 
     ,i,j
     ,evlu1,evlu2,evlu
     ,snr
     ;
     
     Pr:=Generate_IMPolynomial(Pho,x,y):
     N1:=Generate_IMPolynomialmatrix(Pr*P1,1,x,y):
     m:=ImageTools:-Height(N1);
     n:=ImageTools:-Width(N1);
     N11:=Generate_IMPolynomialmatrix(Pr*P1+Nois1,1,x,y):
     evlu:=0:
     evlu1:=0:
     for i from 1 to m do
       for j from 1 to n do
         evlu1:=evlu1+(N1[i,j]-evlu)^2;
       end do;
     end do;
     # the sum of signal
     evlu2:=0:
     for i from 1 to m do
       for j from 1 to n do
         evlu2:=evlu2+(N1[i,j]-N11[i,j])^2;
       end do;
     end do;
     # the sum of noise
     snr:=10*log[10](evlu1/evlu2);
     # the ratio
   
   return snr;
   end proc:

#PP:=InverseFourierTransform(P, inplace = true);
###########################################################################################################################
## In Generate_Polynomial£¬
##     in this algorithm, we get  polynomial FF with the image Create 
##   
## Input:
##     MM---- the image matrix
##
## Output:
##     FF---- the  polynomial we need
###########################################################################################################################

Generate_IMPolynomial:=proc(MM,x,y)
   local
     valu1                                 # the value in computing
     ,m,n                                  # the column dimention and row dimention of the matrix MM
     ,i,j                                  # the parameters in looping
     ,X,Y                                  # the parameters vector of x and y 
     ,FF,MF                                # the polynomial we need with the image matrix
     ;
     n:=ImageTools:-Width(MM);
     m:=ImageTools:-Height(MM);
     X:=Vector(m,i->x^(i-1));
     Y:=Vector(n,i->y^(i-1));
     MF:=Matrix(m,n);
     
     for i from 1 to m do
       for j from 1 to n do
         MF[i,j]:=MM[i,j];
       end do;
     end do;
     FF:=LinearAlgebra:-Transpose(X).MF.Y;
     return expand(FF);
   end proc:

###########################################################################################################################
## In Generate_IMPolynomialmatrix£¬
##     in this algorithm, we get  matrix M with generate image PP be blurred by CC 
##   
## Input:
##     PP---- the polynomial of the real image 
##     CC---- the polynomial for blurring
##
## Output:
##     MM---- the  Create we need
###########################################################################################################################

Generate_IMPolynomialmatrix:=proc(PP,CC,x,y)
   local
     valu1,m,n                              # the value in computing
     ,i,j                                   # the parameters in looping
     ,mm,nn                                 # the column dimention and row dimention of the matrix MM
     ,FF                                    # the polynomial
     ,Vectf_y,Vectf_x                       # the vector
     ,MF                                    # the matrix we need with generate image I
     ; 
    
     FF:=PP*CC;
     m:=degree(FF,x);
     n:=degree(FF,y);
   
     MF:=ImageTools:-Create(m+1,n+1);

     Vectf_y:=CoefficientVector(FF,x);
     
     for i from 1 to m+1 do
       Vectf_x:=CoefficientVector(Vectf_y[i],y);
       mm:=LinearAlgebra:-Dimension(Vectf_x);
       for j from 1 to mm do
         MF[i,j]:=Vectf_x[j];
       end do;
     end do;
         
     return MF;
  end proc:

###########################################################################################################################
## In SD_Generate_IMmatrix£¬
##     in this algorithm, we get  matrix M of image PP  
##   
## Input:
##     PP----  the image 
##
## Output:
##     MF---- the  matrix we need
###########################################################################################################################
SD_Generate_IMmatrix:=proc(PP)
   local
     valu1,m,n                              # the value in computing
     ,i,j                                   # the parameters in looping
     ,MF                                    # the matrix we need of the image 
     ; 
  
     m:=ImageTools:-Height(PP);
     n:=ImageTools:-Width(PP);
     MF:=Matrix(m,n);
     for i from 1 to m do
       for j from 1 to n do
         MF[i,j]:=PP[i,j];
       end do;
     end do;
     return MF;
  end proc:
###########################################################################################################################
## In SD_Generate_matrixIM£¬
##     in this algorithm, we get  image MF of matrix PP  
##   
## Input:
##     PP----  the matrix 
##     evlu---- the coeff we product to PP
##
## Output:
##     MF---- the  Create we need
###########################################################################################################################

SD_Generate_matrixIM:=proc(PP,evlu)
   local
     valu1,m,n                              # the value in computing
     ,i,j                                   # the parameters in looping
     ,MF,MM                                    # the matrix we need of the image 
     ;   
     m:=LinearAlgebra:-RowDimension(PP);
     n:=LinearAlgebra:-ColumnDimension(PP);   
     MF:=ImageTools:-Create(m,n);
     if evlu=1 then
       for i from 1 to m do
         for j from 1 to n do
           MF[i,j]:=PP[i,j];
         end do;
       end do;
     else 
       MM:=evlu*PP;
       for i from 1 to m do
         for j from 1 to n do
           MF[i,j]:=MM[i,j];
         end do;
       end do;
     end if;
         
     return MF;
  end proc:


###########################################################################################################################
## In Generate_matrixZZZZ£¬
##     in this algorithm, we get te Vandermonde matrix with points x1..xk  on the circle of unit and are product by kth unit root: 
##   
## Input:
##     k---- the Row dimension of matrix
##     l---- the Column dimension of matrix
##
## Output:
##     Z---- the Vandermonde matrix
###########################################################################################################################


Generate_matrixZZZZ:=proc(k,l)
   local
     m,n                                 # the number equal to k
     ,i,j                                # the parameters in loops
     ,Vx                                 # the vector of parameters from x1 .. xn
     ,Z                                  # the matrix we need
     ;
     m:=k;
     n:=l;
     Vx:=Vector(m);
     for i from 1 to m do
       Vx[i]:=evalf(cos(2*Pi*(i-1)/m)-I*sin(2*Pi*(i-1)/m));
       # the points
     end do;
     Z:=LinearAlgebra:-VandermondeMatrix(Vx,m,n);
     return Z;
  end proc: