###########Copyright and Notice============== #Copyright since 2001, Bin Han, University of Alberta, #All Rights Reserved #The use of all the following software is permitted for #non-commercial, educational, and research use only. The software #and/or related materials are provided "as-is" without warranty of any #kind including any warranties of performance or merchantability or #fitness for a particular use or purpose or for any purpose #whatsoever, for the licensed product, however used. In no event shall #University of Alberta and/or Bin Han be liable for any damages and/or #costs, including but not limited to incidental or consequential #damages of any kind, including economic damage or injury to property #and lost profits, regardless of whether University of Alberta shall #be advised, have reason to know, or in fact shall know of the #possibility. User bears all risk relating to quality and performance #of the software and/or related materials. #Any use other than non-commercial, educational, or research, or any #redistribution in original or modified form requires prior written #authorization from the copyright holder. #Report errors, mistakes, and bugs to Bin Han at bhan@ualberta.ca. #Send comments/suggestions to Bin Han at bhan@ualberta.ca ##=================================================================== ## MAPLE 6.0 routines for CBC algorithms on constructing 2-dimensional ## biorthogonal multiwavelets with a general 2 by 2 dilation matrix. #This file consists of a set of useful MAPLE subroutines to handle 2D #matrix masks with a general dilation matrix. ### Representation format of matrix mask: # poly(x,y)=sum(Mask[i,j]*x^i*y^j; i, j \in Z): #So poly is a matrix of polynomials with the variables (x,y) reserved #Usage: # the program uses Matrix and Vector from LinearAlgebra. # Do not use matrix and vector from linalg!!! ##D2PrintMask:=proc(poly,choice) ## Print out the matrix mask in various forms ## Input: the matrix of polynomials of the matrix mask. ## Output: Display the matrix mask on screen in various formats. ## and output the matrix mask in array format ## choice<-1, Do not display, output matrix mask in numerical array ## format (6 decimal places) ## choice=-1, Do not display, output matrix mask in exact value array ## choice=0, Display the mask and ouput the exact value form, ## choice=1, Display the mask and ouput the exact value in latex form ## into the file savedlatex ## otherwise, Display and output numerical form (6 Digits) ##D2PrintElement:=proc(poly,xp,yp,choice) ## Print the element of a 2D matrix mask at poistion [xp,yp]. ## Input: the matrix of polynomials of the matrix mask, xp,yp ## Output: Display the element of mask at (xp, yp) on screen. ## choice>=0, exact format, otherwise, numerical to 6 digits ## D2BiorthEqs:=proc(primal_poly,dual_poly,dilation,extrapoly,choice) ## Find the equations of the (bi)orthogonal relation: ## Primalpoly*\bar{Dualpoly} ## is an interpolatory mask w.r.t the lattice Dilation*Z^2. ## Input: primal_poly, dual_poly, dilation, orth, choice ## choice=0, the priaml-dual mask biorthogonal relation ## choice=1, the lowpass--higpass relation completely orthogonal to ## each other ## choice=2, the biorthogonal relation minus given extrapoly ## that is, right side=extrapoly. ## choice>2, compute the relation without the constant ## coefficients. ##D2SumRuleEqs:=proc(poly,dilation,sr,knownvalue,choice) ## Find the equations for sum rules of order sr w.r.t.Dilation*Z^2. ## Input: poly,dilation,sr,choice,knownvalue ## Output: result -> the equations for sum rules of order sr ## choice:=0, then use standard definition of sum rules for ## scalar mask. That is, sum(Mask[i+Mj](i+Mj)^\mu, j\in Z^s) ## is independent of the integer i. ## choice:=1, then each coset equals the knownvalue. That is, ## sum(Mask[i+Mj](M^(-1)i+j)^\mu, j)=knownvalue[\mu]. ## choice=2, the y vector values are known for matrix mask. ## that is, yvect[\mu]=knownvalue[\mu] ## otherwise, we output the equations of sum rules with ## parameterized y vectors (in yvect) ## Do not assume that the sum of the mask is = |det(dilation)| ## The definition of sum rules is defined in Eqs. (1.6) in ## Bin Han, Approximation properties and construction of ## Hermite Interpolants and Biorthogonal Multiwavelets, ## Journal of Approximation Theory, Vol. 110, (2001) 18-53 ##D2CBCAlgorithm:=proc(poly,dilation,sr,choice) ## Given a matrix mask, find the \tilde y vector values by ## the CBC algorithm ## Input: poly,dilation,sr,choice ## Output: result -> the \tilde y vectors in the sum rules ## in Theorem 3.1 in the following paper ## Bin Han, Approximation properties and construction of ## Hermite Interpolants and Biorthogonal Multiwavelets, ## Journal of Approximation Theory, Vol. 110, (2001) 18-53 ## Choice=0, output the moment for scalar mask: ## h^a(\mu):=sum( \tilde a(i+Mj)(M^(-1)i+j)^\mu/\mu!, j) ## Choice=otherwise, output the \tilde y vectors in Theorem 3.1 ##D2MaskOfSym:=proc(sz,dilation,hermiteorder,symgroup,interpolating,cc) ## Parametrize the mask using prescribed symmetry group. ## Input: sz,dilation,hermiteorder,symmetrygroup,interpolating ## Output: a matrix of polynomial for the parameterized matrix mask ## When interpolating=0, it is not an interpolanting mask ## w.r.t. dilation*Z^2, otherwise it is an interpolating matrix mask. ## using the symmetry group. ## cc is the name of the parameters: cc1, cc2, cc3,... ##D2YVectForHermiteMask:=proc(hmord,sr) ##Given hermiteorder and sum rule order, generate the y vector for ##refinable Hermite mask ##D2OutputMask:=proc(poly,sr,choice,filename) ## Output the mask in various form ## choice=0, in numerical format ## otherwise, put the mask in array Mask in parameterize form. ## The parameters are c.number format ##There are also many other useful maple routines in this file. ################################################################# ##Some directly related references: #For symmetry of a refinable function and mask with a dilation matrix # see the papers # Bin Han, Symmetry property and construction of wavelets with a # general dilation matrix, Linear Algebra and Its Applications, # (2001), to appear. # Bin Han, Computing the smoothness exponent of a symmetric # multivariate refinable function, (2001), preprint # Bin Han, Thomas P.-Y. Yu, Bruce Piper, ## Multivariate Refinable Hermite Interpolants, (2002), preprint. #For CBC algorithms on constructing biorthogonal multiwavelets. # CBC stands for Coset By Coset # These program is mainly based on the following paper # Bin Han, Hermite interpolants and biorthogonal multiwavelets # with arbitrary order of vanishing moments (1999) SPIE # Proc. Vol. 3813, pp. 147--161. # The CBC algorithm for the multivariate multiwavelets with # a general dilation matrix was established. # Also, see the following papers for the CBC algorithm # Bin Han, Analysis and Construction of Optimal Multivariate # Biorthogonal Wavelets With Compact Support, SIAM Journal on # Mathematical Analysis, Vol. 31, No.2 (1999/2000), 274--304. # CBC algorithm for interpolatory masks with dilation # matrix 2I_s was first introduced in this paper. # Di-Rong Chen, Bin Han and Sherman D. Riemenschneider, # Construction of Multivariate Biorthogonal Wavelets With # Arbitrary Vanishing Moments, Advances in Computational # Mathematics, Vol. 13 No. 2 (2000), 131-165. # CBC alogrithm for the scalar case with a general dilation # matrix was established. # Bin Han and Rong-Qing Jia, Quincunx Fundamental Refinable # Functions and Quincunx Biorthogonal Wavelets, Mathematics of # Computation, Vol. 71, No. 237, (2002), 165--196. # CBC algorithm for quincunx wavelets was discussed # Bin Han, Construction of multivariate biorthogonal wavelets by # CBC algorithm. Wavelet analysis and multiresolution methods # (Urbana-Champaign, IL, 1999), 105--143, Lecture Notes in Pure # and Appl. Math., 212, Dekker, New York, 2000. # Survey paper on CBC algorithm. # Bin Han, Hermite interpolants and biorthogonal multiwavelets # with arbitrary order of vanishing moments (1999) SPIE # Proc. Vol. 3813, pp. 147--161. # Survey paper on CBC algorithm for 1-dimensional multiwavelets. # Papers can be downloaded at ## http://www.ualberta.ca/~bhan/publ.htm # Program was developed by Bin Han at University of Alberta # Version 1 on January 6, 2002. # Initial Tests of the program have been done on January 6, 2002. # Report bugs, mistakes, errors, comments, suggestions etc. to # Bin Han, bhan@ualberta.ca, htpp://www.ualbert.ca/~bhan #################################################################### ################################################################ ### Main Maple 6.0 routines. ################################################################ read "basicall": default_D2dilation:=Matrix([[2,0],[0,2]]): D2HexSymGroup:=Array(1..12,1..2,1..2, [[[1,0],[0,1]],[[0,-1],[1,-1]],[[-1,1],[-1,0]], [[0,1],[1,0]],[[1,-1],[0,-1]],[[-1,0],[-1,1]], [[-1,0],[0,-1]],[[0,1],[-1,1]],[[1,-1],[1,0]], [[0,-1],[-1,0]],[[-1,1],[0,1]],[[1,0],[1,-1]]], readonly=true): D2FullAxisSymGroup:=Array(1..8,1..2,1..2, [[[1,0],[0,1]],[[-1,0],[0,-1]],[[-1,0],[0,1]],[[1,0],[0,-1]], [[0,1],[1,0]],[[0,-1],[-1,0]],[[0,1],[-1,0]],[[0,-1],[1,0]]], readonly=true): D2XYSymGroup:=Array(1..4,1..2,1..2, [[[1,0],[0,1]],[[-1,0],[0,-1]],[[0,1],[1,0]], [[0,-1],[-1,0]]],readonly=true): D2AxisSymGroup:=Array(1..4,1..2,1..2, [[[1,0],[0,1]],[[-1,0],[0,1]],[[1,0],[0,-1]], [[-1,0],[0,-1]]],readonly=true): D2OrigSymGroup:=Array(1..2,1..2,1..2, [[[1,0],[0,1]],[[-1,0],[0,-1]]],readonly=true): D2NoSymGroup:=Array(1..1,1..2,1..2, [[[1,0],[0,1]]],readonly=true): ##Find the support of a 2D matrix mask. ##Input: the matrix of polynomials of the matrix mask. ##Output: result -> support of mask=[supp[1],supp[2]]*[supp[3],supp[4]] D2Support:=proc(poly) local i,j,multi,POLY,result: result:=Vector(4): multi:=Multiplicity(poly): for i from 1 to multi do for j from 1 to multi do POLY:=expand(poly[i,j]): result[1]:=max(result[1],degree(subs({x=1/x},POLY),x)): result[2]:=max(result[2],degree(POLY,x)): result[3]:=max(result[3],degree(subs({y=1/y},POLY),y)): result[4]:=max(result[4],degree(POLY,y)): od:od: result[1]:=-result[1]: result[3]:=-result[3]: result:=result: end proc: ##Print a 2D matrix mask using various formats. ##Input: the matrix of polynomials of the matrix mask. ##Output: Display the matrix mask on screen in various formats. ## and output the matrix mask in array format ##choice<-1, Do not display, output matrix mask in numerical array ## format (6 decimal places) ##choice=-1, Do not display, output matrix mask in exact value array ##choice=0, Display the mask and ouput the exact value form, ##choice=1, Display the mask and ouput the exact value in latex form ## into the file savedlatex ##otherwise, Display and output numerical form (6 Digits) D2PrintMask:=proc(poly,choice) local i,j,ii,jj,supp,multi,POLY,val,printmask,latexmask,result: multi:=Multiplicity(poly): supp:=D2Support(poly): result:=Array(supp[1]..supp[2],supp[3]..supp[4],1..multi,1..multi): if (choice>=0) then printmask:=Matrix((supp[4]-supp[3]+1)*multi,(supp[2]-supp[1]+1)*multi): fi: for i from 1 to multi do for j from 1 to multi do POLY:=expand(poly[i,j]): for ii from supp[1] to supp[2] do for jj from supp[3] to supp[4] do val:=coeff(coeff(POLY,x,ii),y,jj): if ((choice<(-1))or(choice>1)) then val:=evalf(val, 6): fi: result[ii,jj,i,j]:=val: if (choice>=0) then printmask[(supp[4]-jj)*multi+i,(ii-supp[1])*multi+j]:=val: fi: od:od: od:od: if (choice>=0) then print(printmask, supp): fi: result:=result: if( choice=1) then latexmask:=matrix((supp[4]-supp[3]+1)*multi,(supp[2]-supp[1]+1)*multi,0): for i from 1 to (supp[4]-supp[3]+1)*multi do for j from 1 to (supp[2]-supp[1]+1)*multi do latexmask[i,j]:=printmask[i,j]: od:od: latex(latexmask, savedlatex): fi: printmask:=printmask: end proc: ##Print the element of a 2D matrix mask at poistion [xp,yp]. ##Input: the matrix of polynomials of the matrix mask, xp,yp ##Output: Display the element of the matrix mask on screen. ##choice>=0, exact format, otherwise, numerical to 6 digits D2PrintElement:=proc(poly,xp,yp,choice) local i,j,multi,POLY,val,result: multi:=Multiplicity(poly): result:=Matrix(multi,multi): for i from 1 to multi do for j from 1 to multi do POLY:=expand(poly[i,j]): val:=coeff(coeff(POLY,x,xp),y,yp): if (choice<0) then val:=evalf(val, 6): fi: result[i,j]:=val: od:od: print(result): result:=result: end proc: ##Find the equations of the (bi)orthogonal relation: ## Primalpoly*\bar{Dualpoly} is an interpolatory mask w.r.t ## Dilation*Z^2. ##Input: primal_poly, dual_poly, dilation, orth, choice ##choice=0, the priaml-dual mask biorthogonal relation ##choice=1, the lowpass--higpass relation completely orthogonal to ## each other ##choice=2, the biorthogonal relation minus given extrapoly ## that is, right side=extrapoly. ##choice>2, compute the relation without the constant ## coefficients. D2BiorthEqs:=proc(primal_poly,dual_poly,dilation,extrapoly,choice) local i,j,ii,jj,iii,jjj,m,supp,multi,POLY,result; multi:=Multiplicity(dual_poly): m:=abs(det(dilation)): #absolute value of determinant of dilation matrix POLY:=Matrix(multi, multi): for i from 1 to multi do for j from 1 to multi do POLY[j,i]:=subs({x=1/x, y=1/y},dual_poly[i,j]): od:od: POLY:=.(primal_poly,POLY): if (choice=0) then for i from 1 to multi do POLY[i,i]:=POLY[i,i]-m: od: fi: if (choice=2) then for i from 1 to multi do for j from 1 to multi do POLY[i,j]:=expand(POLY[i,j]-extrapoly[i,j]): od:od: fi: if (choice>2) then for i from 1 to multi do for j from 1 to multi do POLY[i,j]:=expand(POLY[i,j]-coeff(coeff(POLY[i,j],x,0),y,0)): od:od: fi: supp:=D2Support(POLY): result:={}: for i from 1 to multi do for j from 1 to multi do for ii from supp[1] to supp[2] do for jj from supp[3] to supp[4] do iii:=dilation[2,2]*ii-dilation[1,2]*jj: #compute M^(-1)(ii,jj): jjj:=dilation[1,1]*jj-dilation[2,1]*ii: #if (ii,jj) \in dilation*Z^2, then its coefficient must be 0 if((((iii+abs(iii)*m) mod m)=0) and (((jjj+abs(jjj)*m) mod m)=0)) then result:=result union { coeff(coeff(POLY[i,j],x,ii),y,jj) }: fi: od:od:od:od: result:=result: end proc: ##Auxiliary routine for D2Coset D2CosetMap:=proc(dilation) local i,j,ii,jj,kk,m,result: m:=abs(det(dilation)): result:=Array(0..(m*m-1),0): for i from 0 to (m*m-1) do for j from 0 to (m*m-1) do ii:=dilation[2,2]*i-dilation[1,2]*j: #compute M^(-1)(i,j): jj:=dilation[1,1]*j-dilation[2,1]*i: kk:=((ii+abs(ii)*m) mod m)+m*((jj+abs(jj)*m) mod m): result[kk]:=1: od:od: i:=0: for kk from 0 to (m*m-1) do if (result[kk]=1) then result[kk]:=i: i:=i+1: fi: od: result:=result: end proc: ##Given point (i,j), compute the coset of (i,j) w.r.t. dilation*Z^2 D2Coset:=proc(i,j,dilation) local ii,jj,m,result: ii:=dilation[2,2]*i-dilation[1,2]*j: #compute M^(-1)(i,j): jj:=dilation[1,1]*j-dilation[2,1]*i: m:=abs(det(dilation)): result:=((ii+abs(ii)*m) mod m)+m*((jj+abs(jj)*m) mod m): end proc: ##Given matrix E, compute its Kronecker matrix SE of order hmorder ##Input: 2 by 2 matrix E,hmeorder ##Output: result -> (Ex)^\mu/\mu!= \sum( result[\mu,\nu]*x^\nu/\nu!). ## where |\mu|<=hmorder ##choice=0, ouput Array form in 4 indices ##otherwise, output a Matrix with prescribed order. D2KroneckerE:=proc(E,hmorder,choice) local i,j,ii,jj,x,y,val,POLY,result; if (choice=0) then result:=Array(0..hmorder,0..hmorder,0..hmorder,0..hmorder): else result:=Matrix(trunc((hmorder+1)*(hmorder+2)/2)): fi: for i from 0 to hmorder do for j from 0 to i do #[i-j, j] POLY:=expand((E[1,1]*x+E[1,2]*y)^(i-j)*(E[2,1]*x+E[2,2]*y)^j)/((i-j)!*j!): for ii from 0 to hmorder do for jj from 0 to ii do #[ii-jj,jj] val:=coeff(coeff(POLY,x,ii-jj),y,jj)*((ii-jj)!*jj!): if (choice=0) then result[i-j,j,ii-jj,jj]:=val: else result[trunc(i*(i+1)/2)+j+1,trunc(ii*(ii+1)/2)+jj+1]:=val: fi: od:od: od:od: result:=result: end proc: ##Find the moment array of a matrix mask. ##Input: poly,dilation,sr ##Output: result -> the array of the moments of a matrix mask: ##When choice=0, output the whole moments ## Ja[\mu]:=sum(Mask[j]*(M^(-1)j)^\mu, j) ## choice=otherwise, output the coset moments ## Ja[gep, \mu]:=sum(Mask[gep+Mj](M^(-1)gep+j)^\mu, j) D2MomentOfMask:=proc(poly,dilation,sr,choice) local i0,j0,i,j,ii,jj,iii,jjj,m,mm,multi,POLY,supp,val,mask,coset,tmp,result: mm:=det(dilation): m:=abs(mm): multi:=Multiplicity(poly): supp:=D2Support(poly): mask:=Array(supp[1]..supp[2],supp[3]..supp[4],1..multi,1..multi): for i from 1 to multi do for j from 1 to multi do POLY:=expand(poly[i,j]): for ii from supp[1] to supp[2] do for jj from supp[3] to supp[4] do mask[ii,jj,i,j]:=coeff(coeff(POLY,x,ii),y,jj): od:od: od:od: #compute the moment array of a matrix mask if (choice=0) then result:=Array(0..(sr-1),0..(sr-1),1..multi,1..multi): else result:=Array(0..(m-1),0..(sr-1),0..(sr-1),1..multi,1..multi): coset:=D2CosetMap(dilation): fi: for iii from 0 to (sr-1) do for jjj from 0 to iii do #\mu=[iii-jjj,jjj] m:=(iii-jjj)!*(jjj!): #m=\mu! for ii from supp[1] to supp[2] do for jj from supp[3] to supp[4] do i0:=dilation[2,2]*ii-dilation[1,2]*jj: j0:=dilation[1,1]*jj-dilation[2,1]*ii: if(iii=jjj) then val:=1: else val:=(i0/mm)^(iii-jjj): fi: if(jjj<>0) then val:=val*((j0/mm)^jjj): fi: if(choice<>0) then i0:=coset[D2Coset(ii,jj,dilation)]: fi: for i from 1 to multi do for j from 1 to multi do tmp:=val*mask[ii,jj,i,j]/m: if (choice=0) then result[iii-jjj,jjj,i,j]:=result[iii-jjj,jjj,i,j]+tmp: else result[i0,iii-jjj,jjj,i,j]:=result[i0,iii-jjj,jjj,i,j]+tmp: fi: od:od: od:od:od:od: result:=result: end proc: ##Find the equations for sum rules of order sr w.r.t. ## Dilation*Z^2. ##Input: poly,dilation,sr,choice,knownvalue ##Output: result -> the equations for sum rules of order sr ##choice:=0, then use standard definition of sum rules for ## scalar mask. That is, sum(Mask[i+Mj](i+Mj)^\mu, j\in Z^s) ## is independent of the integer i. ##choice:=1, then each coset equals the knownvalue. That is, ## sum(Mask[i+Mj](M^(-1)i+j)^\mu, j)=knownvalue[\mu]. ##choice=2, the y vector values are known for matrix mask. ## that is, yvect[\mu]=knownvalue[\mu] ##otherwise, we output the equations of sum rules with ## parameterized y vectors (in yvect) ##Do not assume that the sum of the mask is = |det(dilation)| ##The definition of sum rules is defined in Eqs. (1.6) in ## Bin Han, Approximation properties and construction of ## Hermite Interpolants and Biorthogonal Multiwavelets, ## Journal of Approximation Theory, Vol. 110, (2001) 18-53 D2SumRuleEqs:=proc(poly,dilation,sr,knownvalue,choice) local i,j,ii,jj,iii,jjj,m,gep,multi,POLY,supp,val,Ja,SE,VAR,tmp,yvect,result; m:=abs(det(dilation)): multi:=Multiplicity(poly): supp:=D2Support(poly): Ja:=D2MomentOfMask(poly,dilation,sr,1): result:=Array(0..sr): for i from 0 to sr do result[i]:={}: od: if ((multi=1)and((choice=0)or(choice=1))) then #multi=1 for ii from 0 to (sr-1) do for jj from 0 to ii do #\mu=[ii-jj,jj] for gep from 0 to (m-1) do if (choice=0) then result[ii+1]:=result[ii+1] union {Ja[gep,ii-jj,jj,1,1]-Ja[0,ii-jj,jj,1,1]}: else result[ii+1]:=result[ii+1] union {Ja[gep,ii-jj,jj,1,1]-knownvalue[ii-jj,jj]}: fi: od:od:od: else SE:=D2KroneckerE(dilation^(-1),sr-1,0): if (choice=2) then yvect:=knownvalue: else yvect:=Array(0..(sr-1),0..(sr-1),1..multi): VAR:={}: jj:=1: for i from 0 to (sr-1) do for j from 0 to i do #[i-j,j] for ii from 1 to multi do yvect[i,j,ii]:=evaln(yvect||jj): VAR:=VAR union { evaln(yvect||jj) }: jj:=jj+1: od:od:od: save VAR, saved: fi: for gep from 0 to (m-1) do for iii from 0 to (sr-1) do for jjj from 0 to iii do #\mu=[iii-jjj,jjj] for i from 1 to multi do tmp:=0: for ii from 0 to iii-jjj do for jj from 0 to jjj do #\nu=[ii,jj] <= \mu for j from 1 to multi do tmp:=tmp+(-1)^(ii+jj)*Ja[gep,ii,jj,j,i]*yvect[iii-jjj-ii,jjj-jj,j]: od:od:od: for ii from 0 to iii do #\nu=[iii-ii,ii], |\nu=|\mu| tmp:=tmp-SE[iii-jjj,jjj,iii-ii,ii]*yvect[iii-ii,ii,i]: od: result[iii+1]:=result[iii+1] union {tmp}: od:od:od:od: fi: for i from 1 to sr do result[0]:=result[0] union convert(result[i],set): od: result:=result: end proc: ##Given a matrix mask, find the \tilde y vector values by ##the CBC algorithm ##Input: poly,dilation,sr,choice ##Output: result -> the \tilde y vectors in the sum rules ## in Theorem 3.1 in the following paper ## Bin Han, Approximation properties and construction of ## Hermite Interpolants and Biorthogonal Multiwavelets, ## Journal of Approximation Theory, Vol. 110, (2001) 18-53 ##Choice=0, output the moment for scalar mask: ## h^a(\mu):=sum( \tilde a(i+Mj)(M^(-1)i+j)^\mu/\mu!, j) ##Choice=otherwise, output the \tilde y vectors in Theorem 3.1 D2CBCAlgorithm:=proc(poly,dilation,sr,choice) local i,ii,iii,j,jj,jjj,k,m,multi,Ja,Ja0MinusmI,SE,tmp,tmp1,VAR,EQ,result: multi:=Multiplicity(poly): m:=abs(det(dilation)): Ja:=D2MomentOfMask(poly,dilation,sr,0): if ((multi=1)and(choice=0)) then if (Ja[0,0,1,1]=0) then printf(The sum of the primal mask must be nonzero!\n): printf(Program is terminated!\n): quit: fi: result:=Array(0..(sr-1),0..(sr-1)): result[0,0]:=1: for ii from 1 to (sr-1) do for jj from 0 to ii do #\mu=[ii-jj,jj] tmp:=0: for iii from 0 to (ii-jj) do for jjj from 0 to jj do #\nu=(iii,jjj) <= \mu if ((iii+jjj)<>ii) then tmp:=tmp+(-1)^(ii-iii-jjj)*result[iii,jjj]*Ja[ii-jj-iii,jj-jjj,1,1]: fi: od:od: result[ii-jj,jj]:=-tmp/Ja[0,0,1,1]: od:od: else result:=Array(0..(sr-1), 0..(sr-1),1..multi): iii:=1: for i from 0 to (sr-1) do for j from 0 to i do for k from 1 to multi do result[i-j,j,k]:=evaln(dualyvector||iii): iii:=iii+1: od:od:od: SE:=D2KroneckerE(dilation^(-1), sr-1,0): Ja0MinusmI:=matrix(multi,multi): for i from 1 to multi do for j from 1 to multi do Ja0MinusmI[i,j]:=Ja[0,0,i,j]: od: Ja0MinusmI[i,i]:=Ja0MinusmI[i,i]-m: od: ##generate y0 vector here tmp1:=kernel(Ja0MinusmI, 'nulldimension'); if (nulldimension <> 1) then printf(There are something wrong with the mask!\n): printf(Program is terminated!\n): quit; fi: tmp:=tmp1[1]: for i from 1 to multi do result[0,0,i]:=tmp[i]: od: #generate all other \tilde y vectors (they must be unique) for iii from 1 to (sr-1) do VAR:={}: EQ:={}: for jjj from 0 to iii do #\mu=(iii-jjj,jjj) for i from 1 to multi do VAR:=VAR union { result[iii-jjj,jjj,i] }: tmp:=-m*result[iii-jjj,jjj,i]: for ii from 0 to (iii-jjj) do for jj from 0 to jjj do #\eta=(ii,jj) <= \mu for j from 0 to (ii+jj) do #\nu=(ii+jj-j,j), |\nu|=|\eta| for k from 1 to multi do tmp:=tmp+SE[ii,jj,ii+jj-j,j]*Ja[iii-jjj-ii,jjj-jj,i,k]*result[ii+jj-j,j,k]: od:od:od:od: EQ:=EQ union {tmp}: od:od: assign(solve(EQ, VAR)): od: fi: result:=result: end proc: ##Parametrize the mask using prescribed symmetry group. ##Input: sz,dilation,hermiteorder,symmetrygroup,interpolating ##Output: a matrix of polynomial for the parameterized matrix mask ##When interpolating=0, it is not an interpolanting mask ##w.r.t. dilation*Z^2, otherwise it is an interpolating matrix mask. ##using the symmetry group ## cc is the name of the parameters: cc1, cc2, cc3,... D2MaskOfSym:=proc(sz,dilation,hermiteorder,symgroup,interpolating,cc) local i,j,ii,jj,iii,jjj,k,num,Eright,Eleft,SEleft,SEright,ind,invdilation,multi,m,VAR,EQ,mask,result: multi:=trunc((hermiteorder+1)*(hermiteorder+2)/2): invdilation:=dilation^(-1): mask:=Array((-sz)..sz,(-sz)..sz,1..multi,1..multi,0): Eleft:=Matrix(2): Eright:=Matrix(2): result:=Matrix(multi): m:=abs(det(dilation)): if (intepolating<>0) then mask[0,0,1..multi,1..multi]:=D2KroneckerE(invdilation,hermiteorder,1): fi: ind:=1: VAR:={}: EQ:={}: ##parameterize the whole mask without symmetry num:=trunc(ArrayNumElems(symgroup)/4): # for ii from 0 to sz do # for jj from 0 to trunc(ii/max(1,(num/4-1))) do #The basic domain for ii from (-sz) to sz do for jj from (-sz) to sz do iii:=dilation[2,2]*ii-dilation[1,2]*jj: jjj:=dilation[1,1]*jj-dilation[2,1]*ii: iii:=iii+abs(iii)*m: jjj:=jjj+abs(jjj)*m: if((interpolating=0)or(((iii mod m)<>0)or((jjj mod m)<>0))) then for i from 1 to multi do for j from 1 to multi do mask[ii,jj,i,j]:=evaln(cc||ind): VAR:= VAR union {evaln(cc||ind)}: ind:=ind+1: od:od: fi: od:od: # ##extend the basic domain to other areas # ##When it folds with the basic domain, set the EQ. num:=trunc(ArrayNumElems(symgroup)/4): for k from 1 to num do #Hexagonal Symmetry Group has 12 elements Eright[1,1]:=symgroup[k][1][1]: Eright[1,2]:=symgroup[k][1][2]: Eright[2,1]:=symgroup[k][2][1]: Eright[2,2]:=symgroup[k][2][2]: Eleft:=.(invdilation,Eright,dilation): SEleft:=D2KroneckerE(Eleft,hermiteorder,1): Eright:=Eright^(-1): SEright:=D2KroneckerE(Eright,hermiteorder,1): # for ii from 0 to sz do # for jj from 0 to trunc(ii/max(1,(num/4-1))) do #The basic domain for ii from (-sz) to sz do for jj from (-sz) to sz do #compute the right side of a(Ek)=S(M^(-1)EM)a(k)S(E^(-1)): for i from 1 to multi do for j from 1 to multi do result[i,j]:=mask[ii,jj,i,j]: od:od: result:=.(SEleft,result,SEright): iii:=symgroup[k,1,1]*ii+symgroup[k,1,2]*jj: jjj:=symgroup[k,2,1]*ii+symgroup[k,2,2]*jj: #if((jjj>=0)and(iii>=0)and(jjj<=trunc(iii/max(1,(num/4-1))))) then if (max(abs(iii),abs(jjj))<=sz) then #[iii,jjj] \in domain for i from 1 to multi do for j from 1 to multi do EQ:=EQ union {result[i,j]-mask[iii,jjj,i,j]}: od:od: else # mask[iii,jjj,1..multi,1..multi]:=result: for i from 1 to multi do for j from 1 to multi do EQ:=EQ union { result[i,j] }: od:od: fi: od:od: od: assign(solve(EQ, VAR)): for i from 1 to multi do for j from 1 to multi do result[i,j]:=0: for ii from (-sz) to sz do for jj from (-sz) to sz do result[i,j]:=result[i,j]+mask[ii,jj,i,j]*(x^ii)*(y^jj): od:od: od:od: result:=result: end proc: ##Given hermiteorder and sum rule order, generate the y vector for ##refinable Hermite mask D2YVectForHermiteMask:=proc(hmord,sr) local i,j,result: if(hmord=0) then result:=Array(0..(sr-1),0..(sr-1),0): result[0,0]:=1: else result:=Array(0..(sr-1),0..(sr-1),1..(trunc((hmord+1)*(hmord+2)/2)),0): for i from 0 to hmord do for j from 0 to i do #\mu=[i-j,j] result[i-j,j,trunc((i+1)*i/2)+j+1]:=1: od:od: fi: result:=result: end proc: ##Output the mask in various form ##choice=0, in numerical format ##otherwise, put the mask in array Mask in parameterize form. ##The parameters are c.number format D2OutputMask:=proc(poly,sr,choice,filename) local i,j,ii,jj,val,multi,supp,Mask: multi:=Multiplicity(poly): supp:=D2Support(poly): interface(quiet=true): writeto(filename): if (choice=0) then #output in numerical format in 16 decimal printf(%d %d %d %d\n,multi,supp[2]-supp[1]+1,supp[4]-supp[3]+1,sr): for ii from supp[1] to supp[2] do for jj from supp[3] to supp[4] do val:=0.0: for i from 1 to multi do for j from 1 to multi do val:=val+abs(evalf(coeff(coeff(poly[i,j],x,ii),y,jj),16)): od:od: if(val>0) then printf(%d %d , ii-supp[1], jj-supp[3]): for i from 1 to multi do for j from 1 to multi do val:=evalf(coeff(coeff(poly[i,j],x,ii),y,jj), 16): printf(%20.16f , val): od:od: printf(\n): fi: od:od: printf(-1\n): else Mask:=array(0..(supp[2]-supp[1]),0..(supp[4]-supp[3]),0..(multi*multi-1)): for ii from supp[1] to supp[2] do for jj from supp[3] to supp[4] do for i from 1 to multi do for j from 1 to multi do val:=evalf(coeff(coeff(poly[i,j],x,ii),y,jj),16): Mask[ii-supp[1],jj-supp[3],(i-1)*multi+j-1]:=val: od:od:od:od: with(codegen, C): C(Mask): fi: writeto(terminal): interface(quiet=false): end proc: #When choice=1, output element in 16 decimal places. #otherwise, output element in exact form. D2OutputTex:=proc(poly,ii,jj,choice) local i,j,multi,Mask: multi:=Multiplicity(poly): Mask:=matrix(multi,multi,0): interface(quiet=true): for i from 1 to multi do for j from 1 to multi do if (choice=1) then Mask[i,j]:=evalf(coeff(coeff(poly[i,j],x,ii),y,jj),16): else Mask[i,j]:=simplify(coeff(coeff(poly[i,j],x,ii),y,jj)): fi: od:od: latex(Mask, OutputTeX): end proc: ##Given a refinement matrix mask poly and a wavelet mask waveletpoly. ##Now the program is used to compute the polynomial matrix ##[psi,psi^T]=sum( int( psi(t)*psi^T(t-(i,j)), t\in R^2)*x^i*y^j). ##So the coefficient of the resulting poly at (0,0) is the L2 norm ##of the wavelet function vector psi. D2BracketPoly:=proc(poly,dilation) local i,ii,iii,j,jj,jjj,m,multi,ind,supp,EQ,VAR,POLY,total,tmp,tmppoly,result; m:=abs(det(dilation)): multi:=Multiplicity(poly): supp:=D2Support(poly): result:=Matrix(multi,multi,0): tmp:=Matrix(multi,multi,0): tmppoly:=Matrix(multi,multi,0): ind:=1: total:=0: VAR:={}: for ii from m*supp[1] to m*supp[2] do for jj from m*supp[3] to m*supp[4] do ##(iii,jjj)=dilation*(ii,jj)^T ? iii:=dilation[1,1]*ii+dilation[1,2]*jj: jjj:=dilation[2,1]*ii+dilation[2,2]*jj: for i from 1 to multi do for j from 1 to multi do result[i,j]:=result[i,j]+evaln(cc||ind)*x^ii*y^jj: tmp[i,j]:=tmp[i,j]+m*evaln(cc||ind)*x^iii*y^jjj: ind:=ind+1: od:od: od:od: for i from 1 to multi do total:=total+coeff(coeff(result[i,i],x,0),y,0): od: EQ:={total-multi}: #The norm of all function=multi tmppoly:=.(poly, result): EQ:=EQ union D2BiorthEqs(tmppoly,poly,dilation,tmp,2): assign(solve(EQ)): for i from 1 to multi do for j from 1 to multi do result[i,j]:=expand(1*result[i,j]): od:od: result:=result: end: D2WavBracketPoly:=proc(wavpoly,dilation,bracketpoly) local i,ii,iii,j,jj,jjj,m,multi,supp,tmp,tmppoly,result; m:=abs(det(dilation)): multi:=Multiplicity(wavpoly): supp:=D2Support(wavpoly): result:=Matrix(multi,multi,0): tmppoly:=Matrix(multi,multi,0): tmp:=Matrix(multi,multi,0): for i from 1 to multi do for j from 1 to multi do tmppoly[i,j]:=subs({x=1/x,y=1/y},wavpoly[j,i]): od:od: tmp:=.(wavpoly, bracketpoly,tmppoly): supp:=D2Support(tmp): m:=det(dilation): for i from 1 to multi do for j from 1 to multi do for ii from supp[1] to supp[2] do for jj from supp[3] to supp[4] do iii:=dilation[2,2]*ii-dilation[1,2]*jj: jjj:=dilation[1,1]*jj-dilation[2,1]*ii: #if (ii,jj) \in dilation*Z^2, if((((iii+abs(iii)*m) mod m)=0)and(((jjj+abs(jjj)*m) mod m)=0)) then result[i,j]:=result[i,j]+coeff(coeff(tmp[i,j],x,ii),y,jj)*x^(iii/m)*y^(jjj/m): fi: od:od:od:od: result:=result: end proc: ##Given the matrix of the polynomials: ##( sum( \int(\psi(t)\psi^*(t-k) dt)*z^k) ##Find the entry norms by take the constants in the diagonal poly. D2WavL2Norm:=proc(wavbracketpoly) local i,multi,result: multi:=Multiplicity(wavbracketpoly): result:=Vector(multi,0): printf("Assume the square sum of all entry norms = %d\n", multi): printf("The entry L2 norms of the function vector = "): for i from 1 to multi do result[i]:=sqrt(coeff(coeff(wavbracketpoly[i,i],x,0),y,0)): printf("%16.12f ", evalf(result[i],24)): od: printf("\n"): result:=result: end proc: ##Convert a poly into a matrix poly D2ConvertToMatrix:=proc(poly) local result: result:=Matrix(1,1,0): result[1,1]:=expand(poly): result:=result: end proc: D2GenPolyXYSym:=proc(sz, num) local i,j,ind,tmp,var,poly: poly:=0: var:={}: tmp:=array(-sz..sz, -sz..sz): for i from -sz to sz do for j from -sz to sz do tmp[i,j]:=0: od:od: ind:=num: for j from 0 to sz do for i from -j to j do tmp[i,j]:=evaln(c||ind): tmp[j,i]:=evaln(c||ind): tmp[-i,-j]:=evaln(c||ind): tmp[-j,-i]:=evaln(c||ind): var := var union { evaln(c||ind) }: ind:=ind+1: od:od: poly:=0: for i from -sz to sz do for j from -sz to sz do poly:=poly+tmp[i,j]*x^i*y^j: od:od: print(var): poly:=expand(poly): end proc: D2GenPoly:=proc(Lx, Hx, Ly, Hy, num) local i,j,ind,poly,var,tmp: poly:=0: var:={}: ind:=num: tmp:=array(Lx..Hx, Ly..Hy): for i from Lx to Hx do for j from Ly to Hy do poly:=poly+evaln(c||ind)*x^i*y^j: var:=var union { evaln(c||ind) }: ind:=ind+1: od:od: print(var): poly:=poly: end proc: D2GenPolyHexSym:=proc(sz, num) local i,j,ind,tmp,poly, var: poly:=0: var:={}: tmp:=array(-sz..sz, -sz..sz): for i from -sz to sz do for j from -sz to sz do tmp[i,j]:=0: od:od: ind:=num: for i from 0 to sz do for j from 0 to trunc(i/2) do tmp[i,j]:=evaln(c||ind): tmp[-i,-j]:=tmp[i,j]: tmp[j,i]:=tmp[i,j]: tmp[-j,-i]:=tmp[i,j]: tmp[i-j,-j]:=tmp[i,j]: tmp[j-i,j]:=tmp[i,j]: tmp[i,i-j]:=tmp[i,j]: tmp[-i,j-i]:=tmp[i,j]: tmp[i-j,i]:=tmp[i,j]: tmp[j-i,-i]:=tmp[i,j]: tmp[j,j-i]:=tmp[i,j]: tmp[-j,i-j]:=tmp[i,j]: var := var union { evaln(c||ind) }: ind:=ind+1: od:od: for i from -sz to sz do for j from -sz to sz do poly:=poly+tmp[i,j]*x^i*y^j: od:od: print(var): poly:=expand(poly): end proc: D2GenPolyFullAxisSym:=proc(sz, num) local i,j,ind,tmp,poly, var: poly:=0: var:={}: tmp:=array(-sz..sz, -sz..sz): for i from -sz to sz do for j from -sz to sz do tmp[i,j]:=0: od:od: ind:=num: for i from 0 to sz do for j from 0 to i do tmp[i,j]:=evaln(c||ind): tmp[i,-j]:=tmp[i,j]: tmp[-i,j]:=tmp[i,j]: tmp[-i,-j]:=tmp[i,j]: tmp[j,i]:=tmp[i,j]: tmp[j,-i]:=tmp[i,j]: tmp[-j,i]:=tmp[i,j]: tmp[-j,-i]:=tmp[i,j]: var := var union { evaln(c||ind) }: ind:=ind+1: od:od: for i from -sz to sz do for j from -sz to sz do poly:=poly+tmp[i,j]*x^i*y^j: od:od: print(var): poly:=expand(poly): end proc: D2GenPolyPointSym:=proc(sx, sy, symx, symy, num) local i,j,ind,tmp,poly, var: poly:=0: var:={}: tmp:=array(symx-sx..sx, symy-sy..sy): for i from (symx-sx) to sx do for j from (symy-sy) to sy do tmp[i,j]:=0: od:od: ind:=num: for i from symx to sx do for j from (symy-sy) to sy do tmp[i,j]:=c(ind): tmp[symx-i,symy-j]:=c(ind): var := var union { c(ind) }: ind:=ind+1: od:od: for i from symx-sx to sx do for j from symy-sy to sy do poly:=poly+tmp[i,j]*x^i*y^j: od:od: print(var): poly:=expand(poly): end proc: