做实习遇到的问题

loveless · 发表于 17.6.2006 16:38:39

<div class='quotetop'>QUOTE(loveless @ 17.06.2006, 16:31 )</div><div class='quotemain'>// Onur G. Guleryuz 1995, 1996, 1997, // University of Illinois at Urbana-Champaign, // Princeton University, // Polytechnic University. #include <stdio.h> #include <stdlib.h> #include "wav_filters_extern.h" #include "wav_gen.h" #include "alloc.h" // Extend data by ofs pixels on both ends using periodic // symmetry. data should be "pointing in" to allow up to data[-ofs] // and data[N-1+ofs]. // ofs>=N is handled. void extend_periodic(float *data,int N,int ofs) { int i,k; k=N-1; for(i=1;i<=ofs;i++) {   data[-i]=data[N-i];   data[k+i]=data[i-1]; } } // Extend data by ofs pixels on both ends using mirror // symmetry. data should be "pointing in" to allow up to data[-ofs] // and data[N-1+ofs]. // // There are two possible "mirrors", i.e., data[-1]=data[1] (mirror on 0) // or data[-1]=data[0] (mirror "between" -1 and 0). phase is either 0, 1, or 2 // to select among these cases. See if(phase==... below. void extend_mirror(float *data,int N,int ofs,int phase) { int i,k; int phase1,phase2; if((phase<0)||(phase>2)) {   printf("extend_mirror: illegal phase\n");   exit(1); } if(phase==2){   // At left boundary the mirror is on 0.   phase1=0;   // At right boundary the mirror is on N-1.   phase2=1; } else {   // phase==0, at left mirror on 0, and at right it's between.   // phase==1, at left mirror between, and at right it's on N-1.   phase1=phase2=phase; } k=N-1;   if(N==1) {   for(i=1;i<=ofs;i++) {   data[-i]=data[0];   data[k+i]=data[0];   } } else   for(i=1;i<=ofs;i++) {   data[-i]=data[i-phase1];   data[k+i]=data[N-phase2-i];   } } // This routine is used in evaluating forward wavelet transforms. // // data of dimension N contains the input. // float *filter contains the Nf filter taps. // ofs is the decimation factor. // // filtering starts at beg and continues until N+beg, yielding N/ofs samples. // data needs to be padded suitably at both ends. // This hoopla about the starting point is in place to ensure // correct boundary processing when evaluating wavelet transforms. // // coef returns the result of the filtering and decimation. // // Viewed as a transform the basis function that generates the // coefficient of index 0 is "<- fliplr(filter) -> (beg)", // i.e., the scalar product obtained by positioning // the flipped filter (formed into a row vector) // with the right side positioned on top of sample beg in data. int filter_n_decimate(float *data,float *coef,int N,float *filter,       int Nf,int beg,int ofs) { int i,j,k; float temp; k=0; for(i=beg;i<N+beg;i+=ofs) {   temp=0;   for(j=i;j>i-Nf;j--)   temp+=data[j]*filter[i-j];   coef[k]=temp;   k++; } return(k); } // This routine is used in evaluating inverse wavelet transforms. // // Upsample and filter to implement a synthesis filter bank. // Basically the inverse of filter_n_decimate() above. // No actual upsampling etc., to avoid zero multiplies. // // coef is the input and contains coefficients resulting from // filtering and decimation by ofs. // // Because this routine is typically called twice for inverse wavelet // transforms, once for the low band and once for high, data is incremented via +=. // Hence data should initially be set to 0. // // Viewed as a transform, filter_n_decimate above starts by putting the // "flipped" forward filter at begf, i.e.,  <- fliplr(ffilter) ->(begf) // The resulting coefficient is at index 0. // Here we are upsampling by ofs and filtering with the inv. filter // time advanced by fbeg so the coefficient at index 0, multiplies the basis function // (-fbeg)<-(ifilter)->(Nf-1-fbeg=beg). // // How is that for confusing? All of this makes sense, // when you write it down carefully. But who has the time ...<!--emo&

--><img src='style_emoticons/<#EMO_DIR#>/smile.gif' border='0' style='vertical-align:middle' alt='smile.gif' /> void upsample_n_filter(float *coef,float *data,int N,float *filter,       int Nf,int beg,int ofs) { int i,j,l,n,p,fbeg; float temp; fbeg=Nf-beg-1; for(i=0;i<N;i++) {   l=0;   n=(i+fbeg)%ofs;   p=(i+fbeg)/ofs;   temp=0;   for(j=n;j<Nf;j+=ofs) {   temp+=coef[p-l]*filter[j];   l++;   }   // += for successive calls for low and high bands.   data+=temp; } } // Used in forward wavelet trf. // // All rows of the image are run through the dual filterbank // given by lp and hp. Results are returned in image "in place". // low freq. coefficients start at 0 and end at N/2, where high freq. // info starts. void filt_n_dec_all_rows(float **image,int Ni,int Nj,float *lp,int Nl,       float *hp,int Nh) { int i,j,ext,ofs; float *data; ext=max(Nl,Nh); data=allocate_1d_float((Nj+2*ext),0); // Point in for required extensions. data+=ext; // offset for the location where high band coefficients // should be copied. ofs=Nj>>1; for(i=0;i<Ni;i++) {     // Copy row.   for(j=0;j<Nj;j++)   data[j]=image[j];   // Extend.   if(PS)   extend_periodic(data,Nj,ext);   else {   // Mirrors at end points.   extend_mirror(data,Nj,ext,2);   }   filter_n_decimate(data,image,Nj,lp,Nl,begflp,2);   filter_n_decimate(data,image+ofs,Nj,hp,Nh,begfhp,2);   } data-=ext; free((void *)data); } // Used in forward wavelet trf. // // Index swapped version of filt_n_dec_all_rows. void filt_n_dec_all_cols(float **image,int Ni,int Nj,float *lp,int Nl,       float *hp,int Nh) { int i,j,ext,ofs; float *data,*data2; ext=max(Nl,Nh); data=allocate_1d_float((Ni+2*ext),0); // Point in for required extensions. data+=ext; // Need second array since cannot access columns directly via pointers. data2=allocate_1d_float(Ni,0); // offset for the location where high band coefficients // should be copied. ofs=Ni>>1; for(j=0;j<Nj;j++) {     // Copy column.   for(i=0;i<Ni;i++)   data=image[j];   // Extend.   if(PS)   extend_periodic(data,Ni,ext);   else {   // Mirrors at end points.   extend_mirror(data,Ni,ext,2);   }   filter_n_decimate(data,data2,Ni,lp,Nl,begflp,2);   filter_n_decimate(data,data2+ofs,Ni,hp,Nh,begfhp,2);   for(i=0;i<Ni;i++)   image[j]=data2;   } data-=ext; free((void *)data); free((void *)data2); } // Used in inverse wavelet trf. // // All rows of the image are run through the dual synthesis filterbank // given by lp and hp having number of taps Nl and Nh respectively. // Results are returned in image "in place". // // Here we actaully need to know the applied shift to the transform, // hence the passed parameters lev and shift_arr. void ups_n_filt_all_rows(float **image,int Ni,int Nj,float *lp,int Nl,       float *hp,int Nh,int lev,int *shift_arr) { int i,j,k,ext,ofs; float *data1,*data2; ext=max(Nl,Nh); ofs=Nj>>1; data1=allocate_1d_float((ofs+2*ext),0); data2=allocate_1d_float((ofs+2*ext),0); data1+=ext;data2+=ext; for(i=0;i<Ni;i++) {   for(j=0;j<ofs;j++) {   k=j+ofs;   // low pass and high pass coefficients.   data1[j]=image[j];image[j]=0;   data2[j]=image[k];image[k]=0;   }   // Take care of the extension at the boundaries.   if(PS) {   extend_periodic(data1,ofs,ext);   extend_periodic(data2,ofs,ext);   }   else {   // shift_arr[lev]=0 OR 1.   // The symmetric banks used are positioned so that on the left side   // the forward lowpass filter has its point of symmetry exactly on top of 0   // and forward high pass filter has its point of symmetry exactly on top of 1.   // This is coordinated by calling filter_n_decimate() with the proper value of beg.   //   // The above is reversed at the right side assuming decimation by 2.   // Hence shift_arr[lev]==0 => at left boundary the mirror is at 0,   // at right boundary it's between, see extend_mirror().   // If shift_arr[lev]==1 then the positions for the low and high are reversed   // since the filters are shifted so if   // shift_arr[lev]==1 => at left boundary the mirror is between,   // at right boundary it's at 0.   extend_mirror(data1,ofs,ext,shift_arr[lev]);   extend_mirror(data2,ofs,ext,1-shift_arr[lev]);   }   // invert low pass.   upsample_n_filter(data1,image,Nj,lp,Nl,begilp,2);   // invert high pass.   upsample_n_filter(data2,image,Nj,hp,Nh,begihp,2); } data1-=ext;data2-=ext; free((void *)data1); free((void *)data2); } // Used in inverse wavelet trf. // // Index swapped version of ups_n_filt_all_rows. void ups_n_filt_all_cols(float **image,int Ni,int Nj,float *lp,int Nl,       float *hp,int Nh,int lev,int *shift_arr) { int i,j,k,ext,ofs; float *data1,*data2,*data3; ext=max(Nl,Nh); ofs=Ni>>1; data1=allocate_1d_float((ofs+2*ext),0); data2=allocate_1d_float((ofs+2*ext),0); data1+=ext;data2+=ext; // Need third array since cannot access columns directly via pointers. data3=allocate_1d_float(Ni,0); for(j=0;j<Nj;j++) {   for(i=0;i<ofs;i++) {   k=i+ofs;   // low pass and high pass coefficients.   data1=image[j];   data2=image[k][j];   data3=data3[k]=0;   }   if(PS) {   extend_periodic(data1,ofs,ext);   extend_periodic(data2,ofs,ext);   }   else {   extend_mirror(data1,ofs,ext,shift_arr[lev]);   extend_mirror(data2,ofs,ext,1-shift_arr[lev]);   }   upsample_n_filter(data1,data3,Ni,lp,Nl,begilp,2);   upsample_n_filter(data2,data3,Ni,hp,Nh,begihp,2);   for(i=0;i<Ni;i++)   image[j]=data3; } data1-=ext;data2-=ext; free((void *)data1); free((void *)data2); free((void *)data3); } [right][snapback]1008971[/snapback][/right] </div> 我其实只需要Wavelet变换的代码。以上的是不是太多了？

loveless · 发表于 17.6.2006 17:09:46

那我该怎么办啊？我把东西都给你，你给我看看行吗？

loveless · 发表于 17.6.2006 17:31:30

好的，我再看下。

loveless · 发表于 17.6.2006 22:22:50

// Onur G. Guleryuz 1995, 1996, 1997, // University of Illinois at Urbana-Champaign, // Princeton University, // Polytechnic University. #include "wav_basic.h" #include "wav_filters_extern.h" #include "alloc.h" #include "wav_gen.h" // Macro that sets filter parameters. // Called via choose_filter(). #define set_filter_param(fname) \ { \ MFLP= fname ## _FLP; \ MFHP= fname ## _FHP; \ MILP= fname ## _ILP; \ MIHP= fname ## _IHP; \ \ /* Filter size. */ \ Nflp= fname ## _Nflp; \ Nfhp= fname ## _Nfhp; \ Nilp= fname ## _Nilp; \ Nihp= fname ## _Nihp; \ \ /* Filter symmetry */ \ PS= fname ## _PS; \ \ /* Filter starting points. */ \ begflp=mfl= fname ## _begflp; \ begfhp=mfh= fname ## _begfhp; \ begilp=mil= fname ## _begilp; \ begihp=mih= fname ## _begihp; \ } // Set wavelet filter parameters depending on name // and tap. See wav_filters.h and wav_filters_extern.h // This routine must be called before doing any transforms // to select which wavelet bank to use. void choose_filter(char name,int tap) { switch (name) { case 'A': switch (tap) { default: set_filter_param(Antonini); } break; case 'H': switch (tap) { default: set_filter_param(Haar); } break; case 'N': switch (tap) { case -1: set_filter_param(Nick_flip); break; default: set_filter_param(Nick); } break; default: switch (tap) { case 9: set_filter_param(D79); break; default: set_filter_param(D8); } } } // 2D "in place" wavelet transform. Input data is in // float **image and output is returned inside image. // // levs>=1 is the number of wavelet levels. // // Ni and Nj are the dimensions of the image. // Ni and Nj should be divisible by 2^levs. // // shift_arr_row and shift_arr_col (each of dimension levs) // contain the amount of shift (0 or 1) // that should be applied to filters in each level for transforms over rows and columns. // The use of shift_arr_row and shift_arr_col is mainly intended for // overcomplete filtering/transforms. // // forw != 0 for a forward transform (inverse is calculated otherwise). // // float *lp of dimension Nl contains the low pass wavelet filter. // float *hp of dimension Nh contains the high pass wavelet filter. // void wav2d_inpl(float **image,int Ni,int Nj,int levs,float *lp,int Nl, float *hp,int Nh,char forw,int *shift_arr_row,int *shift_arr_col) { int i,dimi,dimj; dimi=(forw)?Ni

Ni>>(levs-1)); dimj=(forw)?Nj

Nj>>(levs-1)); for(i=0;i<levs;i++) { if(forw) { // Do the rows. // Adjust filter parameters for overcomplete filtering. begflp=mfl+shift_arr_row; begfhp=mfh-shift_arr_row; filt_n_dec_all_rows(image,dimi,dimj,lp,Nl,hp,Nh); // Now the columns. begflp=mfl+shift_arr_col; begfhp=mfh-shift_arr_col; filt_n_dec_all_cols(image,dimi,dimj,lp,Nl,hp,Nh); } else { // Rows. begilp=mil+shift_arr_row[levs-1-i]; begihp=mih-shift_arr_row[levs-1-i]; ups_n_filt_all_rows(image,dimi,dimj,lp,Nl,hp,Nh,levs-1-i,shift_arr_row); // Columns. begilp=mil+shift_arr_col[levs-1-i]; begihp=mih-shift_arr_col[levs-1-i]; ups_n_filt_all_cols(image,dimi,dimj,lp,Nl,hp,Nh,levs-1-i,shift_arr_col); } dimi=(forw)?(dimi>>1)dimi<<1); dimj=(forw)?(dimj>>1)dimj<<1); } } // 2-D "in place" wavelet packet transform. // See wav2d_inpl() for an explanation of parameters. // // This routine will grow a full packet tree for levels upto LL_ONLY_AFTER_LEV // and just a wavelet tree afterward, i.e., if levs>LL_ONLY_AFTER_LEV // we grow wavelets only over LL band after LL_ONLY_AFTER_LEV. // LL_ONLY_AFTER_LEV is set in wav_gen.h Just set it to some large number // if you want packets all the way. // // Ni and Nj should be divisible by 2^levs. // shift_arr_row and shift_arr_col contain the amount of shift (0 or 1) // that should be applied to filters in each level. // This is mainly intended for overcomplete filtering. // See wav2d_inpl() comments. // void wavpack2d_inpl(float **image,int Ni,int Nj,int levs,float *lp,int Nl, float *hp,int Nh,char forw,int *shift_arr_row,int *shift_arr_col) { int i,k,l,dimi,dimj; int p,q,rcnt,ccnt; float **buffer; buffer=allocate_2d_float(Ni,Nj,0); dimi=(forw)?NiNi>>(levs-1)); dimj=(forw)?Nj:(Nj>>(levs-1)); for(i=0;i<levs;i++) { // Number of row bands. rcnt=Ni/dimi; // Number of column bands. ccnt=Nj/dimj; if(forw&&(i>=LL_ONLY_AFTER_LEV)) { // Grow wavelets only over LL band after LL_ONLY_AFTER_LEV. rcnt=ccnt=1; } else if((!forw)&&(levs-1-i>=LL_ONLY_AFTER_LEV)) { // Inverse wavelets only over LL band after LL_ONLY_AFTER_LEV. rcnt=ccnt=1; } // Do the next level of the wavelet tree over all of // the bands we are supposed to grow on. for(p=0;p<rcnt;p++) for(q=0;q<ccnt;q++) { for(k=0;k<dimi;k++) for(l=0;l<dimj;l++) buffer[k][l]=image[p*dimi+k][q*dimj+l]; if(forw) { // Do the rows. // Adjust filter parameters for overcomplete filtering. begflp=mfl+shift_arr_row; begfhp=mfh-shift_arr_row; filt_n_dec_all_rows(buffer,dimi,dimj,lp,Nl,hp,Nh); // Now the columns. begflp=mfl+shift_arr_col; begfhp=mfh-shift_arr_col; filt_n_dec_all_cols(buffer,dimi,dimj,lp,Nl,hp,Nh); } else { // Rows. begilp=mil+shift_arr_row[levs-1-i]; begihp=mih-shift_arr_row[levs-1-i]; ups_n_filt_all_rows(buffer,dimi,dimj,lp,Nl,hp,Nh,levs-1-i,shift_arr_row); // Columns. begilp=mil+shift_arr_col[levs-1-i]; begihp=mih-shift_arr_col[levs-1-i]; ups_n_filt_all_cols(buffer,dimi,dimj,lp,Nl,hp,Nh,levs-1-i,shift_arr_col); } for(k=0;k<dimi;k++) for(l=0;l<dimj;l++) image[p*dimi+k][q*dimj+l]=buffer[k][l]; } dimi=(forw)?(dimi>>1):(dimi<<1); dimj=(forw)?(dimj>>1):(dimj<<1); } free_2d_float(buffer,Ni); } // Forward complex wavelet transform using Kingsbury's paper. // Please check with Nick Kingsbury's paper to confirm that // I am doing things correctly here. // // Input is in float **im, and the results are returned in 4 // **float's, trf[0],...,trf[3]. // // Ni, Nj multiples of 2^levs. // levs>=1. // See wav2d_inpl() for an explanation of parameters. // void complex_wav_forw(float **im,float ***trf,int Ni,int Nj,int levs) { int i,j,k,l; int Nl,Nh,dimi,dimj,hdimi,hdimj; float *lp,*hp; int shift_arr_r[MAX_ARR_SIZE],shift_arr_c[MAX_ARR_SIZE]; float tmp; // Choose Antonini. choose_filter('A',0); lp=MFLP;Nl=Nflp; hp=MFHP;Nh=Nfhp; // First level, fully overcomplete resulting 4x redundancy. // Results are in trf. for(i=0;i<4;i++) { for(k=0;k<Ni;k++) for(l=0;l<Nj;l++) trf[k][l]=im[k][l]; // Set shifts as specified by Kingsbury's paper. if(i/2==0) shift_arr_r[0]=1; else shift_arr_r[0]=0; if(i%2==0) shift_arr_c[0]=1; else shift_arr_c[0]=0; // Transform. wav2d_inpl(trf,Ni,Nj,1,lp,Nl,hp,Nh,1,shift_arr_r,shift_arr_c); } // Now grow complex wavelets over the LL band // using Kingsbury's orthogonal bank. if(levs>1) { for(i=0;i<4;i++) { dimi=Ni/2; dimj=Nj/2; for(j=1;j<levs;j++) { if(i/2==0) { // Regular orth. bank. choose_filter('N',0); } else { // Flipped orth. bank (see Kingsbury). // Flipped is same as inverse. choose_filter('N',-1); } lp=MFLP;Nl=Nflp; hp=MFHP;Nh=Nfhp; // Transform over rows. filt_n_dec_all_rows(trf,dimi,dimj,lp,Nl,hp,Nh); if(i%2==0) { // Regular. choose_filter('N',0); } else { // Flipped. choose_filter('N',-1); } lp=MFLP;Nl=Nflp; hp=MFHP;Nh=Nfhp; // Transform over columns. filt_n_dec_all_cols(trf,dimi,dimj,lp,Nl,hp,Nh); dimi>>=1; dimj>>=1; } } } // Generate the complex wavelet coefficients. hdimi=Ni/(1<<levs); hdimj=Nj/(1<<levs); for(k=0;k<Ni;k++) for(l=0;l<Nj;l++) { // After this computation the complex wavelet coefficients // are given by trf[0]+jtrf[2] and trf[3]+jtrf[1]; tmp=trf[0][k][l]+trf[3][k][l]; trf[3][k][l]=trf[0][k][l]-trf[3][k][l]; trf[0][k][l]=tmp; tmp=trf[1][k][l]+trf[2][k][l]; // It is possible that the below is 2 - 1. trf[2][k][l]=trf[1][k][l]-trf[2][k][l]; trf[1][k][l]=tmp; } } // Inverse complex wavelet transform. // trf is modified in intermediate computations. // Ni, Nj multiples of 2^levs // levs>=1. // // See complex_wav_forw() comments. // float **complex_wav_inv(float ***trf,int Ni,int Nj,int levs) { int i,j,k,l; int Nl,Nh,dimi,dimj,hdimi,hdimj; float *lp,*hp; int shift_arr_r[MAX_ARR_SIZE],shift_arr_c[MAX_ARR_SIZE]; float tmp; float **im; // Clean the filter shift arrays. for(i=0;i<1<<levs;i++) shift_arr_r=shift_arr_c=0; // Get back the intermediate wavelet coefficients from the // complex wavelet coefficients. hdimi=Ni/(1<<levs); hdimj=Nj/(1<<levs); for(k=0;k<Ni;k++) for(l=0;l<Nj;l++) { // Prior to this computation the complex wavelet coefficients // are given by trf[0]+jtrf[2] and trf[3]+jtrf[1]; tmp=(trf[0][k][l]+trf[3][k][l])/2; trf[3][k][l]=(trf[0][k][l]-trf[3][k][l])/2; trf[0][k][l]=tmp; tmp=(trf[1][k][l]+trf[2][k][l])/2; trf[2][k][l]=(trf[1][k][l]-trf[2][k][l])/2; trf[1][k][l]=tmp; } if(levs>1) { for(i=0;i<4;i++) { dimi=Ni>>(levs-1); dimj=Nj>>(levs-1); for(j=1;j<levs;j++) { if(i/2==0) { // Regular inverse, it just happens that orth. inverses are flipped. choose_filter('N',0); } else { // Flipped inverse. choose_filter('N',-1); } lp=MILP;Nl=Nilp; hp=MIHP;Nh=Nihp; // Inverse transform over rows. ups_n_filt_all_rows(trf,dimi,dimj,lp,Nl,hp,Nh,levs-1-j,shift_arr_r); if(i%2==0) { // Regular inverse. choose_filter('N',0); } else { // Flipped inverse. choose_filter('N',-1); } lp=MILP;Nl=Nilp; hp=MIHP;Nh=Nihp; // Inverse transform over columns. ups_n_filt_all_cols(trf,dimi,dimj,lp,Nl,hp,Nh,levs-1-j,shift_arr_c); dimi<<=1; dimj<<=1; } } } // Initialize final result with zeros. im=allocate_2d_float(Ni,Nj,1); // Choose Antonini. choose_filter('A',0); lp=MILP;Nl=Nilp; hp=MIHP;Nh=Nihp; // Now the first level. for(i=0;i<4;i++) { // Set shifts as specified by Kingsbury's paper. if(i/2==0) shift_arr_r[0]=1; else shift_arr_r[0]=0; if(i%2==0) shift_arr_c[0]=1; else shift_arr_c[0]=0; // Inverse transform. wav2d_inpl(trf,Ni,Nj,1,lp,Nl,hp,Nh,0,shift_arr_r,shift_arr_c); // Accumulate the results of 4x redundancy. for(k=0;k<Ni;k++) for(l=0;l<Nj;l++) im[k][l]+=trf[k][l]; } // Normalize results by 4 to account for 4x redundancy accumulation. for(k=0;k<Ni;k++) for(l=0;l<Nj;l++) im[k][l]/=4; return(im); } // Forward complex wavelet packet transform. // Ni, Nj multiples of 2^levs. // levs>=1. // // See complex_wav_forw() comments. // See wavpack2d_inpl() comments for LL_ONLY_AFTER_LEV. // // Not fully tested. Please check with literature. // void complex_wav_pack_forw(float **im,float ***trf,int Ni,int Nj,int levs) { int i,j,k,l; int Nl,Nh,dimi,dimj,hdimi,hdimj; float *lp,*hp; int shift_arr_r[MAX_ARR_SIZE],shift_arr_c[MAX_ARR_SIZE]; float tmp,**buffer; int p,q,rcnt,ccnt; // Choose Antonini. choose_filter('A',0); lp=MFLP;Nl=Nflp; hp=MFHP;Nh=Nfhp; // First level, fully overcomplete resulting 4x redundancy. // Results are in trf. for(i=0;i<4;i++) { for(k=0;k<Ni;k++) for(l=0;l<Nj;l++) trf[k][l]=im[k][l]; // Set shifts as specified by Kingsbury's paper. if(i/2==0) shift_arr_r[0]=1; else shift_arr_r[0]=0; if(i%2==0) shift_arr_c[0]=1; else shift_arr_c[0]=0; // Transform. wav2d_inpl(trf,Ni,Nj,1,lp,Nl,hp,Nh,1,shift_arr_r,shift_arr_c); } // Now grow using Kingsbury's orthogonal bank. if(levs>1) { buffer=allocate_2d_float(Ni/2,Nj/2,0); for(i=0;i<4;i++) { dimi=Ni/2; dimj=Nj/2; for(j=1;j<levs;j++) { // Number of row/column bands. rcnt=Ni/dimi; ccnt=Nj/dimj; if(j>=LL_ONLY_AFTER_LEV) { // Grow wavelets only over LL band after LL_ONLY_AFTER_LEV. rcnt=ccnt=1; } for(p=0;p<rcnt;p++) for(q=0;q<ccnt;q++) { // Copy data to buffer. for(k=0;k<dimi;k++) for(l=0;l<dimj;l++) buffer[k][l]=trf[p*dimi+k][q*dimj+l]; if(i/2==0) { // Regular orth. bank. choose_filter('N',0); } else { // Flipped orth. bank (see Kingsbury). // Flipped is same as inverse. choose_filter('N',-1); } lp=MFLP;Nl=Nflp; hp=MFHP;Nh=Nfhp; // Transform over rows. filt_n_dec_all_rows(buffer,dimi,dimj,lp,Nl,hp,Nh); if(i%2==0) { // Regular. choose_filter('N',0); } else { // Flipped. choose_filter('N',-1); } lp=MFLP;Nl=Nflp; hp=MFHP;Nh=Nfhp; // Transform over columns. filt_n_dec_all_cols(buffer,dimi,dimj,lp,Nl,hp,Nh); // Copy results to trf. for(k=0;k<dimi;k++) for(l=0;l<dimj;l++) trf[p*dimi+k][q*dimj+l]=buffer[k][l]; } dimi>>=1; dimj>>=1; } } free_2d_float(buffer,Ni/2); } // Generate the complex wavelet coefficients. hdimi=Ni/(1<<levs); hdimj=Nj/(1<<levs); for(k=0;k<Ni;k++) for(l=0;l<Nj;l++) { // After this computation the complex wavelet coefficients // are given by trf[0]+jtrf[2] and trf[3]+jtrf[1]; tmp=trf[0][k][l]+trf[3][k][l]; trf[3][k][l]=trf[0][k][l]-trf[3][k][l]; trf[0][k][l]=tmp; tmp=trf[1][k][l]+trf[2][k][l]; trf[2][k][l]=trf[1][k][l]-trf[2][k][l]; trf[1][k][l]=tmp; } } // Inverse complex wavelet packet transform. // trf is modified in intermediate computations. // Ni, Nj multiples of 2^levs // levs>=1. // // See complex_wav_pack_forw() comments. // float **complex_wav_pack_inv(float ***trf,int Ni,int Nj,int levs) { int i,j,k,l; int Nl,Nh,dimi,dimj,hdimi,hdimj; float *lp,*hp; int shift_arr_r[MAX_ARR_SIZE],shift_arr_c[MAX_ARR_SIZE]; float tmp; float **im,**buffer; int p,q,rcnt,ccnt; // Clean the filter shift arrays. for(i=0;i<1<<levs;i++) shift_arr_r=shift_arr_c=0; // Get back the intermediate wavelet coefficients from the // complex wavelet coefficients. hdimi=Ni/(1<<levs); hdimj=Nj/(1<<levs); for(k=0;k<Ni;k++) for(l=0;l<Nj;l++) { // Prior to this computation the complex wavelet coefficients // are given by trf[0]+jtrf[2] and trf[3]+jtrf[1]; tmp=(trf[0][k][l]+trf[3][k][l])/2; trf[3][k][l]=(trf[0][k][l]-trf[3][k][l])/2; trf[0][k][l]=tmp; tmp=(trf[1][k][l]+trf[2][k][l])/2; trf[2][k][l]=(trf[1][k][l]-trf[2][k][l])/2; trf[1][k][l]=tmp; } if(levs>1) { buffer=allocate_2d_float(Ni/2,Nj/2,0); for(i=0;i<4;i++) { dimi=Ni>>(levs-1); dimj=Nj>>(levs-1); for(j=1;j<levs;j++) { // Number of row/column bands. rcnt=Ni/dimi; ccnt=Nj/dimj; if(j>=LL_ONLY_AFTER_LEV) { // Grow wavelets only over LL band after LL_ONLY_AFTER_LEV. rcnt=ccnt=1; } for(p=0;p<rcnt;p++) for(q=0;q<ccnt;q++) { // Copy data to buffer. for(k=0;k<dimi;k++) for(l=0;l<dimj;l++) buffer[k][l]=trf[p*dimi+k][q*dimj+l]; if(i/2==0) { // Regular inverse, it just happens that orth. inverses are flipped. choose_filter('N',0); } else { // Flipped inverse. choose_filter('N',-1); } lp=MILP;Nl=Nilp; hp=MIHP;Nh=Nihp; // Inverse transform over rows. ups_n_filt_all_rows(buffer,dimi,dimj,lp,Nl,hp,Nh,levs-1-j,shift_arr_r); if(i%2==0) { // Regular inverse. choose_filter('N',0); } else { // Flipped inverse. choose_filter('N',-1); } lp=MILP;Nl=Nilp; hp=MIHP;Nh=Nihp; // Inverse transform over columns. ups_n_filt_all_cols(buffer,dimi,dimj,lp,Nl,hp,Nh,levs-1-j,shift_arr_c); // Copy results to trf. for(k=0;k<dimi;k++) for(l=0;l<dimj;l++) trf[p*dimi+k][q*dimj+l]=buffer[k][l]; } dimi<<=1; dimj<<=1; } } free_2d_float(buffer,Ni/2); } // Initialize final result with zeros. im=allocate_2d_float(Ni,Nj,1); // Choose Antonini. choose_filter('A',0); lp=MILP;Nl=Nilp; hp=MIHP;Nh=Nihp; // Now the first level. for(i=0;i<4;i++) { // Set shifts as specified by Kingsbury's paper. if(i/2==0) shift_arr_r[0]=1; else shift_arr_r[0]=0; if(i%2==0) shift_arr_c[0]=1; else shift_arr_c[0]=0; // Inverse transform. wav2d_inpl(trf,Ni,Nj,1,lp,Nl,hp,Nh,0,shift_arr_r,shift_arr_c); // Accumulate the results of 4x redundancy. for(k=0;k<Ni;k++) for(l=0;l<Nj;l++) im[k][l]+=trf[k][l]; } // Normalize results by 4 to account for 4x redundancy accumulation. for(k=0;k<Ni;k++) for(l=0;l<Nj;l++) im[k][l]/=4; return(im); }

loveless · 发表于 17.6.2006 22:25:13

我实在看得不行了，这个是不是正宗的Wavelet Transformation? 那个很多的，是不是很多多余的？

德国开元华人社区开元周游

新闻看板

找房租房

招聘求职

经验分享

德国生活

买卖市场

留德新生

留学德国

开元交友

心情水吧

美食天地

分类广告

做实习遇到的问题

请您关注我们在:

站点信息

站点统计| 举报| Archiver| 手机版| 小黑屋

Powered by Discuz! X3.2 © 2001-2014 Comsenz Inc.

GMT+1, 13.11.2024 22:22

关于我们|Apps

() 开元网