|
|
// Onur G. Guleryuz 1995, 1996, 1997,<br />// University of Illinois at Urbana-Champaign,<br />// Princeton University,<br />// Polytechnic University.<br /><br />#include <stdio.h><br />#include <stdlib.h><br />#include "wav_filters_extern.h"<br />#include "wav_gen.h"<br />#include "alloc.h"<br /><br />// Extend data by ofs pixels on both ends using periodic <br />// symmetry. data should be "pointing in" to allow up to data[-ofs]<br />// and data[N-1+ofs]. <br />// ofs>=N is handled.<br /><br />void extend_periodic(float *data,int N,int ofs)<br /><br />{<br /> int i,k;<br /><br /> k=N-1;<br /> for(i=1;i<=ofs;i++) {<br /><br /> data[-i]=data[N-i];<br /> data[k+i]=data[i-1];<br /> }<br />}<br /><br />// Extend data by ofs pixels on both ends using mirror <br />// symmetry. data should be "pointing in" to allow up to data[-ofs]<br />// and data[N-1+ofs]. <br />//<br />// There are two possible "mirrors", i.e., data[-1]=data[1] (mirror on 0) <br />// or data[-1]=data[0] (mirror "between" -1 and 0). phase is either 0, 1, or 2<br />// to select among these cases. See if(phase==... below. <br /> <br />void extend_mirror(float *data,int N,int ofs,int phase)<br /><br />{<br /> int i,k;<br /> int phase1,phase2;<br /><br /> if((phase<0)||(phase>2)) {<br /><br /> printf("extend_mirror: illegal phase\n");<br /> exit(1);<br /> }<br /><br /> if(phase==2){<br /> // At left boundary the mirror is on 0.<br /> phase1=0;<br /><br /> // At right boundary the mirror is on N-1.<br /> phase2=1;<br /> }<br /> else {<br /><br /> // phase==0, at left mirror on 0, and at right it's between.<br /> // phase==1, at left mirror between, and at right it's on N-1.<br /> phase1=phase2=phase;<br /> }<br /><br /> k=N-1;<br /> <br /> if(N==1) {<br /><br /> for(i=1;i<=ofs;i++) {<br /><br /> data[-i]=data[0];<br /> data[k+i]=data[0];<br /> }<br /> }<br /> else<br /> for(i=1;i<=ofs;i++) {<br /><br /> data[-i]=data[i-phase1];<br /> data[k+i]=data[N-phase2-i];<br /> }<br />}<br /><br />// This routine is used in evaluating forward wavelet transforms.<br />//<br />// data of dimension N contains the input.<br />// float *filter contains the Nf filter taps.<br />// ofs is the decimation factor.<br />//<br />// filtering starts at beg and continues until N+beg, yielding N/ofs samples. <br />// data needs to be padded suitably at both ends.<br />// This hoopla about the starting point is in place to ensure <br />// correct boundary processing when evaluating wavelet transforms.<br />//<br />// coef returns the result of the filtering and decimation.<br />// <br />// Viewed as a transform the basis function that generates the<br />// coefficient of index 0 is "<- fliplr(filter) -> (beg)", <br />// i.e., the scalar product obtained by positioning <br />// the flipped filter (formed into a row vector) <br />// with the right side positioned on top of sample beg in data.<br /><br />int filter_n_decimate(float *data,float *coef,int N,float *filter,<br /> int Nf,int beg,int ofs)<br /><br />{<br /> int i,j,k;<br /> float temp;<br /><br /> k=0;<br /> for(i=beg;i<N+beg;i+=ofs) {<br /><br /> temp=0;<br /> for(j=i;j>i-Nf;j--)<br /> temp+=data[j]*filter[i-j];<br /><br /> coef[k]=temp;<br /> k++;<br /> }<br /><br /> return(k);<br />}<br /><br />// This routine is used in evaluating inverse wavelet transforms.<br />//<br />// Upsample and filter to implement a synthesis filter bank.<br />// Basically the inverse of filter_n_decimate() above.<br />// No actual upsampling etc., to avoid zero multiplies.<br />//<br />// coef is the input and contains coefficients resulting from <br />// filtering and decimation by ofs. <br />//<br />// Because this routine is typically called twice for inverse wavelet<br />// transforms, once for the low band and once for high, data is incremented via +=.<br />// Hence data should initially be set to 0.<br />//<br />// Viewed as a transform, filter_n_decimate above starts by putting the <br />// "flipped" forward filter at begf, i.e., <- fliplr(ffilter) ->(begf)<br />// The resulting coefficient is at index 0.<br />// Here we are upsampling by ofs and filtering with the inv. filter<br />// time advanced by fbeg so the coefficient at index 0, multiplies the basis function<br />// (-fbeg)<-(ifilter)->(Nf-1-fbeg=beg).<br />// <br />// How is that for confusing? All of this makes sense,<br />// 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' /><!--endemo--><br /><br />void upsample_n_filter(float *coef,float *data,int N,float *filter,<br /> int Nf,int beg,int ofs)<br /><br />{<br /> int i,j,l,n,p,fbeg;<br /> float temp;<br /><br /> fbeg=Nf-beg-1;<br /><br /> for(i=0;i<N;i++) {<br /><br /> l=0;<br /> n=(i+fbeg)%ofs;<br /> p=(i+fbeg)/ofs;<br /><br /> temp=0;<br /> for(j=n;j<Nf;j+=ofs) {<br /><br /> temp+=coef[p-l]*filter[j];<br /> l++;<br /> }<br /><br /> // += for successive calls for low and high bands. <br /> data+=temp;<br /> }<br />}<br /><br />// Used in forward wavelet trf.<br />//<br />// All rows of the image are run through the dual filterbank<br />// given by lp and hp. Results are returned in image "in place".<br />// low freq. coefficients start at 0 and end at N/2, where high freq.<br />// info starts.<br />void filt_n_dec_all_rows(float **image,int Ni,int Nj,float *lp,int Nl,<br /> float *hp,int Nh)<br /><br />{<br /> int i,j,ext,ofs;<br /> float *data;<br /><br /> ext=max(Nl,Nh);<br /> data=allocate_1d_float((Nj+2*ext),0);<br /> // Point in for required extensions.<br /> data+=ext;<br /><br /> // offset for the location where high band coefficients<br /> // should be copied.<br /> ofs=Nj>>1;<br /><br /> for(i=0;i<Ni;i++) {<br /> <br /> // Copy row.<br /> for(j=0;j<Nj;j++)<br /> data[j]=image[j];<br /><br /> // Extend.<br /> if(PS)<br /> extend_periodic(data,Nj,ext);<br /> else {<br /> // Mirrors at end points.<br /> extend_mirror(data,Nj,ext,2);<br /> }<br /><br /> filter_n_decimate(data,image,Nj,lp,Nl,begflp,2);<br /> filter_n_decimate(data,image+ofs,Nj,hp,Nh,begfhp,2);<br /> }<br /><br /> data-=ext;<br /> free((void *)data);<br /> }<br /><br />// Used in forward wavelet trf.<br />//<br />// Index swapped version of filt_n_dec_all_rows.<br />void filt_n_dec_all_cols(float **image,int Ni,int Nj,float *lp,int Nl,<br /> float *hp,int Nh)<br /><br />{<br /> int i,j,ext,ofs;<br /> float *data,*data2;<br /><br /> ext=max(Nl,Nh);<br /> data=allocate_1d_float((Ni+2*ext),0);<br /> // Point in for required extensions.<br /> data+=ext;<br /><br /> // Need second array since cannot access columns directly via pointers.<br /> data2=allocate_1d_float(Ni,0);<br /><br /> // offset for the location where high band coefficients<br /> // should be copied.<br /> ofs=Ni>>1;<br /><br /> for(j=0;j<Nj;j++) {<br /> <br /> // Copy column.<br /> for(i=0;i<Ni;i++)<br /> data=image[j];<br /><br /> // Extend.<br /> if(PS)<br /> extend_periodic(data,Ni,ext);<br /> else {<br /> // Mirrors at end points.<br /> extend_mirror(data,Ni,ext,2);<br /> }<br /><br /> filter_n_decimate(data,data2,Ni,lp,Nl,begflp,2);<br /> filter_n_decimate(data,data2+ofs,Ni,hp,Nh,begfhp,2);<br /><br /> for(i=0;i<Ni;i++)<br /> image[j]=data2;<br /> }<br /><br /> data-=ext;<br /> free((void *)data);<br /> free((void *)data2);<br /> }<br /><br />// Used in inverse wavelet trf.<br />//<br />// All rows of the image are run through the dual synthesis filterbank<br />// given by lp and hp having number of taps Nl and Nh respectively. <br />// Results are returned in image "in place".<br />//<br />// Here we actaully need to know the applied shift to the transform,<br />// hence the passed parameters lev and shift_arr.<br /><br />void ups_n_filt_all_rows(float **image,int Ni,int Nj,float *lp,int Nl,<br /> float *hp,int Nh,int lev,int *shift_arr)<br /><br />{<br /> int i,j,k,ext,ofs;<br /> float *data1,*data2;<br /><br /> ext=max(Nl,Nh);<br /> ofs=Nj>>1;<br /><br /> data1=allocate_1d_float((ofs+2*ext),0);<br /> data2=allocate_1d_float((ofs+2*ext),0);<br /> data1+=ext;data2+=ext;<br /><br /> for(i=0;i<Ni;i++) { <br /><br /> for(j=0;j<ofs;j++) {<br /><br /> k=j+ofs;<br /> // low pass and high pass coefficients.<br /> data1[j]=image[j];image[j]=0;<br /> data2[j]=image[k];image[k]=0;<br /> }<br /><br /> // Take care of the extension at the boundaries.<br /> if(PS) {<br /><br /> extend_periodic(data1,ofs,ext); <br /> extend_periodic(data2,ofs,ext); <br /> }<br /> else {<br /><br /> // shift_arr[lev]=0 OR 1.<br /> // The symmetric banks used are positioned so that on the left side<br /> // the forward lowpass filter has its point of symmetry exactly on top of 0<br /> // and forward high pass filter has its point of symmetry exactly on top of 1. <br /> // This is coordinated by calling filter_n_decimate() with the proper value of beg.<br /> //<br /> // The above is reversed at the right side assuming decimation by 2. <br /> // Hence shift_arr[lev]==0 => at left boundary the mirror is at 0, <br /> // at right boundary it's between, see extend_mirror().<br /><br /> // If shift_arr[lev]==1 then the positions for the low and high are reversed<br /> // since the filters are shifted so if<br /> // shift_arr[lev]==1 => at left boundary the mirror is between, <br /> // at right boundary it's at 0.<br /> extend_mirror(data1,ofs,ext,shift_arr[lev]); <br /> extend_mirror(data2,ofs,ext,1-shift_arr[lev]); <br /> }<br /><br /> // invert low pass.<br /> upsample_n_filter(data1,image,Nj,lp,Nl,begilp,2);<br /> // invert high pass.<br /> upsample_n_filter(data2,image,Nj,hp,Nh,begihp,2);<br /> }<br /><br /> data1-=ext;data2-=ext;<br /> free((void *)data1);<br /> free((void *)data2);<br />}<br /><br />// Used in inverse wavelet trf.<br />//<br />// Index swapped version of ups_n_filt_all_rows.<br />void ups_n_filt_all_cols(float **image,int Ni,int Nj,float *lp,int Nl,<br /> float *hp,int Nh,int lev,int *shift_arr)<br /><br />{<br /> int i,j,k,ext,ofs;<br /> float *data1,*data2,*data3;<br /><br /> ext=max(Nl,Nh);<br /> ofs=Ni>>1;<br /><br /> data1=allocate_1d_float((ofs+2*ext),0);<br /> data2=allocate_1d_float((ofs+2*ext),0);<br /> data1+=ext;data2+=ext;<br /><br /> // Need third array since cannot access columns directly via pointers.<br /> data3=allocate_1d_float(Ni,0);<br /><br /> for(j=0;j<Nj;j++) { <br /><br /> for(i=0;i<ofs;i++) {<br /><br /> k=i+ofs;<br /> // low pass and high pass coefficients.<br /> data1=image[j];<br /> data2=image[k][j];<br /> data3=data3[k]=0;<br /> }<br /><br /> if(PS) {<br /><br /> extend_periodic(data1,ofs,ext); <br /> extend_periodic(data2,ofs,ext); <br /> }<br /> else {<br /><br /> extend_mirror(data1,ofs,ext,shift_arr[lev]); <br /> extend_mirror(data2,ofs,ext,1-shift_arr[lev]); <br /> }<br /><br /> upsample_n_filter(data1,data3,Ni,lp,Nl,begilp,2);<br /> upsample_n_filter(data2,data3,Ni,hp,Nh,begihp,2);<br /><br /> for(i=0;i<Ni;i++)<br /> image[j]=data3;<br /> }<br /><br /> data1-=ext;data2-=ext;<br /> free((void *)data1);<br /> free((void *)data2);<br /> free((void *)data3);<br />}<br /> |
|
楼主: loveless
狗仔卡
发表于 17.6.2006 08:43:40
显身卡
楼主
rogram Files\MATLAB704\toolbox\matlab\general\mex.m is a case-insensitive match and will be used instead. You can improve the performance of your code by using exact name matches and we therefore recommend that you update your usage accordingly. Alternatively, you can disable this warning using warning('off','MATLAB:dispatcher:InexactMatch').<br /> <br />Select a compiler: <br />[1] Lcc C version 2.4 in C:\