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



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


##########################################################################################################################
## Firstly,we must get the code of approximate cofactor u(x)utilzing Bezout matrix 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 cofactor u(x,y) 
## directly using svd decomposition.
## At last we can get the code of completely.
##########################################################################################################################


########################################################################################################################### 
## In BZ_Get_1Dcofactors,
##     in this algorithm, we may compute the cofactor of f (degree(f)>=degree(g))with the approximate 1D_GCD degree utilizing 
##      leading sub Bezout matrix: 
##                       
## Input:
##     f,g ----the 1-D polynomials
##      deg----the 1D_gcd's degree
##       
## Output 
##     u-----the cofactor of f
##########################################################################################################################


BZ_Get_1Dcofactors:=proc(f,g,x,deg)
   # deg is the degree of approximate gcd
   local 
     m,n                 # the degree of f and g
     ,MM,NN              # the subsylvester matrix
     ,V                  # the singularvalues of the sylvestermatrix
     ,W,W1,VV               # the singular vecters
     ,i,j
     ,V1,V2
     ,u,v                # the polynomials we need
     ,st
     ;
     st:=time():
     m:=degree(f,x);
     n:=degree(g,x);
     NN:=SD_Get_subBezoutmatrix(f,g,deg,x);
     VV:=LinearAlgebra:-SingularValues(NN,output=['Vt']);
     W1:=LinearAlgebra:-HermitianTranspose(LinearAlgebra:-Row(VV,m-deg+1));
     st:=time()-st:
     W:=Get_LPBezoutmatrix(f,deg,x).W1; 
     u:=0;
     # computing the u where f/u is the gcd
     for i from 1 to m-deg+1 do
       u:=u+W(i)*x^(i-1);
     end do;
     u:=expand((u/lcoeff(u))/norm(expand(u/lcoeff(u)),2));
     userinfo(1,S_stls,print(`the cofactors' coeff time is`),print('st'=st)); 
     return u;
   end proc:


########################################################################################################################### 
## In BZ_Get_gcddegree,
##     in this algorithm, we may compute the degree of the 1D degree utilizing leading subBezout matrix: 
##                       
## Input:
##     f,g ----the 1-D polynomials
##       
## Output 
##     deg-----the degre of the gcd of f and g
##########################################################################################################################

BZ_ClGet_gcddegree:=proc(f,g,x,e)
   # e is the tolerance times
   local 
     m,n,k                 # the degree of f and g
     ,MM,NN                # the subsylvester matrix
     ,_S,_SS               # the singularvalues of the sylvestermatrix
     ,i
     ,deg                  # the dgree we need
     ,st
     ,ff,gg
     ;
     st:=time():
     m:=degree(f,x);
     n:=degree(g,x);
     ff:=f/norm(expand(f),2):
     gg:=g/norm(expand(g),2):
     deg:=min(m,n);
     deg:=deg-8;
     # the deg must smaller than min(m,n) to be less
     MM:=SD_Get_subBezoutmatrix(f,g,deg,x);
     _S:=LinearAlgebra:-SingularValues(LinearAlgebra:-Transpose(MM),output='list');
     k:=m-deg+1;
     while _S[k]<10^(-e) and deg<min(m,n) do
       k:=k-1;       
     end do:
     deg:=m-k:
     st:=time()-st:
     userinfo(1,S_stls,print(`Getting degree using time is: `),print('st1'=st),print('degy'=_S)); 
     return deg;
   end proc:

########################################################################################################################### 
## In SD_Get_subBezoutmatrix,
##     in this algorithm, we may compute the leading subBezout matrix of f and g with given gcd degree: 
##                       
## Input:
##     f,g ----the 1-D polynomials
##     deg ----the degree
##       
## Output 
##     SM-----the subBezoutmatrix we need
##########################################################################################################################

SD_Get_subBezoutmatrix:=proc(f,g,deg,x)
   # e is the tolerance times
   local 
     m,n              # the degree of f and g
     ,BM              # the return
     ,V1,V2           # the coefficient vector of f and g 
     ,i,j,k,s         # the parameters for loops
     ,evlu            # the value for computing
     ,st
     ;
     st:=time():
     m:=degree(f,x);
     n:=degree(g,x);
     V1:=CoefficientVector(f, x);
     V2:=CoefficientVector(g, x):
     BM:=Matrix(m-deg+1,m-deg+1):
     for i from 1 to m-deg+1 do
       for j from 1 to i do
	 k:=m-i;
	 # k=m-i as B[j,i]=b(m-j+1,m-i+1)
	 s:=0;
	 # s=0..k
	 evlu:=0;
	 while m-j+s+1<m+1 and s <k+1 do          # <=m means not muiti the zero
	   if m-j+s+1>n then                      # >n means only have one terms or both zero
	     if k-s<n+1 then
	       evlu:=evlu+V1[m-j+s+1+1]*V2[k-s+1];
	     else
	       evlu:=evlu+0;
	     end if;
	   else
	     if k-s<n+1 then                      # <n means every one can if not means only one
	       evlu:=evlu+V1[m-j+s+1+1]*V2[k-s+1]-V1[k-s+1]*V2[m-j+s+1+1];
	     else
	       evlu:=evlu-V1[k-s+1]*V2[m-j+s+1+1];
	     end if;
	   end if;
	   s:=s+1;
	 end do;
	 BM[j,i]:=evlu;
       end do;
     end do; 
     for i from 2 to m-deg+1 do
       for j from 1 to i-1 do
         BM[i,j]:=BM[j,i]:
       end do;
     end do; 
     st:=time()-st:
     userinfo(3,S_stls,print(`the cofactors' coeff time is`),print('st'=st));
     return BM;
   end proc:

########################################################################################################################### 
## In Get_LPBezoutmatrix,
##     in this algorithm, we may compute the leading principal leaind sub Bezout matrix of f and 1 with given gcd degree: 
##                       
## Input:
##     f ----the 1-D polynomials
##     deg ----the degree
##       
## Output 
##     SM-----the LPBezoutmatrix we need
##########################################################################################################################


Get_LPBezoutmatrix:=proc(f,deg,x)
   # e is the tolerance times
   local 
     m,n              # the degree of f and g
     ,BM              # the return
     ,V1,V2           # the coefficient vector of f and g 
     ,i,j,k,s         # the parameters for loops
     ,evlu            # the value for computing
     ,st
     ;
     st:=time():
     m:=degree(f,x);
     n:=degree(g,x);
     V1:=CoefficientVector(f, x);
     BM:=Matrix(m-deg+1,m-deg+1):
     for i from 1 to m-deg+1 do
       for j from i to m-deg+1 do
         BM[i,j]:=V1[m+1-j+i];
       end do;
     end do;     
     st:=time()-st:
     userinfo(3,S_stls,print(`the cofactors' coeff time is`),print('st'=st));
     return BM;
   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 cofactor 1th
###########################################################################################################################
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
     ,m1,n1                                         # the degree of f with x and y
     ,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();

     
     st:= time();
     m:=LinearAlgebra:-RowDimension(M_f):
     n:=LinearAlgebra:-ColumnDimension(M_f):
     m1:=m-1;
     n1:=n-1;     
     X :=Vector(m, i->x^(i-1)); 
     Y :=Vector(n, i->x^(i-1));
     st:=time()-st:
     userinfo(1,S_stls,print(`the begin AB time is`),print('st'=st));
     st:=time():
     Z3  :=Generate_matrixZZZZ(1,m);
     ZZ3 :=Generate_matrixZZZZ(1,n);
     F1 := Row(Z3, 1).M_f.Y;
     F2 := Row(Z3, 1).M_g.Y;
     F1:=F1/norm(expand(F1),infinity):F2:=F2/norm(expand(F2),infinity):
     degy:=BZ_ClGet_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),infinity):FF2:=FF2/norm(expand(FF2),infinity):
     degx:=BZ_ClGet_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));
     Z  := Generate_matrixZZZZ(m1+1-degx,m);
     ZZ := Generate_matrixZZZZ(n1+1-degy,n);
     A:=Matrix(m1+1-degx,n1+1-degy);
     B:=Matrix(m1+1-degx,n1+1-degy);
     st:= time()-st;
     userinfo(3,S_stls,print(`the interpolation vector  is `),print('st1'=Z,ZZ));
     st:= time();
     deg:=degy;
     for i from 1 to m1+1-degx do
       st2:= time();
       F1 := Row(Z, i).M_f.Y;
       F2 := Row(Z, i).M_g.Y;
       st2:= time()-st2;
       userinfo(3,S_stls,print(`the computing time 1D f and g  is `),print('st1'=st2));
       st2:= time();
       F:=BZ_Get_1Dcofactors(F1,F2,x,deg);
       st2:= time()-st2;
       userinfo(3,S_stls,print(`the cofactor computing time 1D f and g  is `),print('st1'=st2));
       st2:= time();
       if F=1 then
         for j from 1 to n1+1-degy do
           A[i,j] := 1;
         end do;
       else
         for j from 1 to n1+1-degy do
	   A[i,j] := eval(F,x=evalf(cos(2*Pi*(j-1)/(n1+1-degy))-I*sin(2*Pi*(j-1)/(n1+1-degy))));
         end do;
       end if;
       st2:= time()-st2;
       userinfo(3,S_stls,print(`thecofactor computing time 1D f and g  is `),print('st1'=st2));
     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 cofactor's computing time is `),print('st1'=st));
     st:= time();
     deg:=degx;
     for i from 1 to n1+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)));
       st2:= time()-st2;
       userinfo(1,S_stls,print(`the evaluation time 1D f and g  is `),print('st1'=st2));
       st2:= time();
       FF:=BZ_Get_1Dcofactors(FF1,FF2,x,deg);
       st2:= time()-st2;
       userinfo(3,S_stls,print(`the cofactors computing time 1D f and g  is `),print('st1'=st2));
       st2:= time();
       if FF=1 then
         for j from 1 to m1+1-degx do
           B[j,i] := 1;
         end do;
       else
	 for j from 1 to m1+1-degx do
	   B[j,i] := eval(FF,x=evalf(cos(2*Pi*(j-1)/(m1+1-degx))-I*sin(2*Pi*(j-1)/(m1+1-degx))));
         end do;
       end if;
       st2:= time()-st2;
       userinfo(3,S_stls,print(`the cofactor1's computing time 1D f and g  is `),print('st1'=st2),print('deg'=degx));
       userinfo(3,S_stls,print(`generate the 1D_GCD of f and g :`),print('PP'=FF));
     end do; 
     st:= time()-st;
     userinfo(1,S_stls,print(`the second time cofactor'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_cofactor1,
##     in this algorithm, get the co-prime part of the first argument : 
##
## Input:
##     M1,M2---coeff matrices of 2D polynomials
##     e------ tolerance
##     
##       
## Output 
##     CC ----the co-prime part of the first argument
###########################################################################################################################
Get_cofactor1:=proc(M1,M2, e :: numeric)
   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 Inverse fast FourierTransform 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);
     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; 
     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 Ap2D_gcd,
##     in this algorithm, we may compute the coeff matrix of 2-D GCD of M1 and M2 and coeff matrix of cofactor u(x,y) : 
##
## Input:
##     M1,M2----the 2-D polynomials' coeff matrix
##     e    ---- the tolerance
##     u   ---- the matrix of coeff of u(x,y)( zero or not)
##       
## Output 
##     P3,CC-----the coeff matrix of 2-D polynomials we need
###########################################################################################################################
Ap2D_gcd:=proc( M1,M2,e,u)
   local 
     P1,P2,P3,P4                         # the matrixes of P and InverseFourierTransform its as defined in the paper 
     ,Q                                  # the polynomials of approximate 2D-gcd of f and g      
     ,m1,n1,m2,n2,m,n,mm,nn              # the row and column of Matrix A&B defined in the paper
     ,deg                                # the degree of GCD
     ,evlu                               # the parameters for computing
     ,i,j                                # the parameters in computing for loop statement
     ,st,st1                             # the time beginning
     ,CC,DD 
     ,xdeg,ydeg
     :
     st:= time():
     m1:=LinearAlgebra:-RowDimension(M1):
     m2:=LinearAlgebra:-RowDimension(M2):
     n1:=LinearAlgebra:-ColumnDimension(M1):
     n2:=LinearAlgebra:-ColumnDimension(M2):
     if u=0 then
       st1:= time(): 
       CC:=Get_cofactor1(M1,M2,e);
       st1:= time()-st1;
       userinfo(1,S_stls,print(`we get co-fator 1 using computing time is `),print('st1'=st1));
     else
       CC:=u;
     end if;
     st1:= time():     
     P3:=SD_Dividen(M1,CC);
     st1:= time()-st1;
     userinfo(1,S_stls,print(`we get GCD by dividen using computing time is `),print('st1'=st1));
     st:= time()-st;
     userinfo(1,S_stls,print(`we get gcd using computing time is `),print('st1'=st));
     return P3,CC;
end proc:  
########################################################################################################################### 
## In RGB_ImageDeconvolution,
##   we will deconvolution the origin version from the input two RGB distorted versions
##
##   Input  Pho1,Pho2---- the blurred image versions(RGB )
##              e    ---- the accuracy
##
##   Output  Pho    ---- the version we deconvolution
########################################################################################################################### 
RGB_ImageDeconvolution:=proc(Pho1,Pho2,e)
   local
     m1,n1,m2,n2,m3,n3                          # the Hight and width of these two version
     ,img_r1,img_g1,img_b1                      # the input Pho1's R G B
     ,img_r2,img_g2,img_b2                      # the input Pho2's R G B
     ,m,n                                       # the output dimension
     ,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
     ,aevlu,aevlu1,aevlu2,aevlu3                # the value for computing
     ,Pho                                       # the output we need
     ,st
     ;
     st:=time():
     m1:=ImageTools:-Height(Pho1);
     n1:=ImageTools:-Width(Pho1);
     m2:=ImageTools:-Height(Pho2);
     n2:=ImageTools:-Width(Pho2);
     # the dimension of these two picture
     img_r1:= ImageTools:-Create(m1, n1, (i,j) -> Pho1[i,j,1]):
     img_g1:= ImageTools:-Create(m1, n1, (i,j) -> Pho1[i,j,2]):
     img_b1:= ImageTools:-Create(m1, n1, (i,j) -> Pho1[i,j,3]):
     img_r2:= ImageTools:-Create(m2, n2, (i,j) -> Pho2[i,j,1]):
     img_g2:= ImageTools:-Create(m2, n2, (i,j) -> Pho2[i,j,2]):
     img_b2:= ImageTools:-Create(m2, n2, (i,j) -> Pho2[i,j,3]):
     # and there R G B to grey photo
     st:=time()-st:
     userinfo(1,S_stls,print(`Before we compute gcd  using computing time is `),print('st1'=st));
     # get one of gcd
     img_r,u1:=SD_ImageDeconvolution(img_r1,img_r2,e,0);
     # get one of gcd
     img_g,u2:=SD_ImageDeconvolution(img_g1,img_g2,e,0);
     # get one of gcd
     img_b,u3:=SD_ImageDeconvolution(img_b1,img_b2,e,0);
     Pho:= ImageTools:-CombineLayers(img_r, img_g, img_b):
     return Pho,u1;
  end proc:


########################################################################################################################### 
## In SD_ImageDeconvolution,
##   we will deconvolution the origin version from the input three distorted versions
##
##   Input  Pho1,Pho2---- the blurred image versions( Grey)
##              e    ---- the accuracy
##
##   Output  Pho    ---- the version we deconvolution
########################################################################################################################### 
SD_ImageDeconvolution:=proc(Pho1,Pho2,e,uu)
   local
     m1,n1,m2,n2,m3,n3                          # the Hight and width of these three version
     ,m,n                                       # the output dimension
     ,f,g,h                                     # the 2D polynomials
     ,p,u                                       # the polynomial is 2D gcd of f and g and cofactor
     ,M1,M2 
     ,evlu                                      # the value for computing
     ,aevlu                                     # the value for computing
     ,Pho                                       # the output we need
     ,st
     ;
     st:=time():
     m1:=ImageTools:-Height(Pho1):
     n1:=ImageTools:-Width(Pho1):
    # the dimension of these two picture
     M1:=SD_Generate_IMmatrix(Pho1):
     M2:=SD_Generate_IMmatrix(Pho2):
     st:=time()-st:
     userinfo(1,S_stls,print(`Before we compute gcd  using computing time is `),print('st1'=st));
     p,u:=Ap2D_gcd(M1,M2,e,0,uu);
     # get one of gcd
     st:=time():
     evlu:=Pho1[m1,n1]:
     Pho:=SD_Generate_matrixIM(p,(6)*evlu):
     st:=time()-st:
     userinfo(1,S_stls,print(`After we compute gcd  using computing time is `),print('st1'=st));
     return Pho,u;
  end proc:
########################################################################################################################### 
## In SSD_ImageDeconvolution,
##   we will deconvolution the origin version from the input three distorted versions
##
##   Input  Pho1,Pho2,Pho3---- the blurred image versions(RGB or Grey)
##              e    ---- the accuracy
##
##   Output  Pho    ---- the version we deconvolution
########################################################################################################################### 
SSD_ImageDeconvolution:=proc(Pho1,Pho2,Pho3,e)
   local
     m1,n1                                       # the Hight and width of these three version
     ,f,g,h                                     # the 2D polynomials
     ,p,u                                       # the polynomial is 2D gcd of f and g and cofactor
     ,M1,M2,M3
     ,evlu,evlu1,evlu2,evlu3                    # the value for computing
     ,aevlu,aevlu1,aevlu2,aevlu3                # the value for computing
     ,Pho                                       # the output we need
     ,st
     ;
     st:=time():
     M1:=SD_Generate_IMmatrix(Pho1):
     M2:=SD_Generate_IMmatrix(Pho2):
     M3:=SD_Generate_IMmatrix(Pho3):
     evlu1:= 7/10;#RandomTools[Generate](integer(range=1..10))/10;#1/2;#
     evlu2:= 3/5;#RandomTools[Generate](integer(range=1..10))/10;#1/5;#
     evlu3:= 1;#RandomTools[Generate](integer(range=1..10))/10;#1/5;#
     aevlu1:=1/5;# RandomTools[Generate](integer(range=1..10))/10;#1;#
     aevlu2:= 3/5;#RandomTools[Generate](integer(range=1..10))/10;#1;#
     aevlu3:=2/5;# RandomTools[Generate](integer(range=1..10))/10;#7/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,u:=Ap2D_gcd(evlu1*M1+evlu2*M3+evlu3*M3,aevlu1*M1+aevlu2*M2+aevlu3*M3,e,0);
     # get the gcd
     st:=time():
     m1:=ImageTools:-Height(Pho1);
     n1:=ImageTools:-Width(Pho1);
     evlu:=Pho1[m1,n1]:
     Pho:=SD_Generate_matrixIM(p,6*evlu):
     st:=time()-st:
     userinfo(1,S_stls,print(`After we compute gcd  using computing time is `),print('st1'=st));
     return Pho,u;
  end proc:

########################################################################################################################### 
## In SSD_Blurred,
##   we will blur the image Pho into two version from P1 and P2 ,P3
##
##   Input  Pho --- the image need to be blurred
##          P1,P2,P3--- the blur polynomial
##
##   Output  Pho1,Pho2 ----the versions we need
########################################################################################################################### 
SRGB_Blurred:=proc(Pho,Pr1,Pr2,Pg1,Pg2,Pb1,Pb2,Noisr1,Noisr2,Noisg1,Noisg2,Noisb1,Noisb2)
   local
     m,n                                        # the dimension Pho 
     ,img_r,img_g,img_b                         # the trip-pho pictures
     ,Nr1,Nr2,Nr3,Ng1,Ng2,Ng3,Nb1,Nb2,Nb3       # the two versions of R G B 
     ,Pho1,Pho2                                 # the return
     ;
     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
     Nr1:=SD_Blurred(img_r,Pr1,Noisr1):
     Nr2:=SD_Blurred(img_r,Pr2,Noisr2):
     #return img_r,Nr1,Nr2;
     # the blurred R versions
     Ng1:=SD_Blurred(img_g,Pg1,Noisg1):
     Ng2:=SD_Blurred(img_g,Pg2,Noisg2):
     #return img_g,Ng1,Ng2;
     # the blurred G versions
     Nb1:=SD_Blurred(img_b,Pb1,Noisb1):
     Nb2:=SD_Blurred(img_b,Pb2,Noisb2):
     #return img_b,Nb1,Nb2;
     # the blurred B versions
     Pho1:= ImageTools:-CombineLayers(Nr1, Ng1, Nb1):
     Pho2:= ImageTools:-CombineLayers(Nr2, Ng2, Nb2):
     # the combine RGB versions
   return Pho1,Pho2;
   end proc:
########################################################################################################################### 
## In SD_Blurred,
##   we will blur the image Pho into two version from P1 and P2 ,P3
##
##   Input  Pho --- the image need to be blurred
##          P1--- the blur polynomial
##          noise1 ---- the noise
##
##   Output  Pho1 ----the versions we need
########################################################################################################################### 
SD_Blurred:=proc(Pho,P1,Nois1)
   local
     PP                                         # the polynomials of these trip-pho
     ,Pho1                                      # the return
     ;    
     # the grey versions of Pho
     PP:=Generate_IMPolynomial(Pho,x,y):
     Pho1:=Generate_IMPolynomialmatrix(PP*P1+Nois1,1,x,y):     
     # the blurred  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
########################################################################################################################### 
SNR:=proc(Pho,Pr1,Pr2,Pg1,Pg2,Pb1,Pb2,Noisr1,Noisr2,Noisg1,Noisg2,Noisb1,Noisb2)
   local
     m,n                           # the dimension Pho 
     ,img_r,img_g,img_b             # the trip-pho pictures
     ,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,Pr1,Noisr1);
     snr2r:=SD_SNR(img_r,Pr2,Noisr2);
     snr1g:=SD_SNR(img_g,Pg1,Noisg1);
     snr2g:=SD_SNR(img_g,Pg2,Noisg2);
     snr1b:=SD_SNR(img_b,Pb1,Noisb1);
     snr2b:=SD_SNR(img_b,Pb2,Noisb2);
     # 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:
########################################################################################################################### 
## 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 FFT 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 FFT 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:


#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: