wav_trf.c

matlab+dct+数字水印

			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[i],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[i],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[i],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[i][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[i].
	for(i=0;i<4;i++) {

		for(k=0;k<Ni;k++)
			for(l=0;l<Nj;l++)
				trf[i][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[i],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[i][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[i].
						for(k=0;k<dimi;k++)
							for(l=0;l<dimj;l++)
								trf[i][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[i]=shift_arr_c[i]=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[i][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[i].
						for(k=0;k<dimi;k++)
							for(l=0;l<dimj;l++)
								trf[i][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[i],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[i][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);
}