DIP程序.rar_8-bit Histogram_dct idct_fft ifft_图像处理_数字图像处理 灰度

/*############################################################################# * 文件名：dct.c * 功能： 离散余弦变换 * modified by PRTsinghua@hotmail.com #############################################################################*/ #include "wm.h" #include "dct.h" #define INVROOT2 0.7071067814 #define SWAP(A, B) {double t = A; A = B; B = t;} int N; int M; double *dct_NxN_tmp = NULL; double *dct_NxN_costable = NULL; int dct_NxN_log2N = 0; static const unsigned int JPEG_lumin_quant_table[NJPEG][NJPEG] = { {16, 11, 10, 16, 24, 40, 51, 61}, {12, 12, 14, 19, 26, 58, 60, 55}, {14, 13, 16, 24, 40, 57, 69, 56}, {14, 17, 22, 29, 51, 87, 80, 62}, {18, 22, 37, 56, 68, 109, 103, 77}, {24, 35, 55, 64, 81, 104, 113, 92}, {49, 64, 78, 87, 103, 121, 120, 101}, {72, 92, 95, 98, 112, 100, 103, 99}}; static const unsigned int JPEG_chromin_quant_table[NJPEG][NJPEG] = { {17, 18, 24, 47, 99, 99, 99, 99}, {18, 21, 26, 66, 99, 99, 99, 99}, {24, 26, 56, 99, 99, 99, 99, 99}, {47, 66, 99, 99, 99, 99, 99, 99}, {99, 99, 99, 99, 99, 99, 99, 99}, {99, 99, 99, 99, 99, 99, 99, 99}, {99, 99, 99, 99, 99, 99, 99, 99}, {99, 99, 99, 99, 99, 99, 99, 99}}; // 初始化余弦数组 static void initcosarray() { int i,group,base,item,nitems,halfN; double factor; dct_NxN_log2N = -1; do{ dct_NxN_log2N++; if ((1<<dct_NxN_log2N)>N){ fprintf(stderr, "dct_NxN: %d not a power of 2\n", N); exit(1); } }while((1<<dct_NxN_log2N)<N); if (dct_NxN_costable) free(dct_NxN_costable); dct_NxN_costable = malloc(N * sizeof(double)); #ifdef DEBUG if(!dct_NxN_costable){ fprintf(stderr, "Unable to allocate C array\n"); exit(1); } #endif halfN=N/2; for(i=0;i<=halfN-1;i++) dct_NxN_costable[halfN+i]=4*i+1; for(group=1;group<=dct_NxN_log2N-1;group++){ base= 1<<(group-1); nitems=base; factor = 1.0*(1<<(dct_NxN_log2N-group)); for(item=1; item<=nitems;item++) dct_NxN_costable[base+item-1]=factor*dct_NxN_costable[halfN+item-1]; } for(i=1;i<=N-1;i++) dct_NxN_costable[i] = 1.0/(2.0*cos(dct_NxN_costable[i]*M_PI/(2.0*N))); } // 初始化N×N的dct void init_dct_NxN(int width, int height) { #ifdef DEBUG if (width != height || width <= 0) { fprintf(stderr, "init_dct_NxN(): dimensions out of range\n"); exit(1); } #endif if (dct_NxN_tmp && M != height) free(dct_NxN_tmp); N = width; M = height; dct_NxN_tmp = malloc(height * sizeof(double)); #ifdef DEBUG if (!dct_NxN_tmp) { fprintf(stderr, "init_dct_NxN(): failed to allocate memory\n"); exit(1); } #endif initcosarray(); } // 位反转 static void bitrev(double *f, int len) { int i,j,m; if (len<=2) return; /* No action necessary if n=1 or n=2 */ j=1; for(i=1; i<=len; i++){ if(i<j) SWAP(f[j-1], f[i-1]); m = len>>1; while(j>m){ j=j-m; m=(m+1)>>1; } j=j+m; } } // 和反转 static void inv_sums(double *f) { int stepsize,stage,curptr,nthreads,thread,step,nsteps; for(stage=1; stage <=dct_NxN_log2N-1; stage++){ nthreads = 1<<(stage-1); stepsize = nthreads<<1; nsteps = (1<<(dct_NxN_log2N-stage)) - 1; for(thread=1; thread<=nthreads; thread++){ curptr=N-thread; for(step=1; step<=nsteps; step++){ f[curptr] += f[curptr-stepsize]; curptr -= stepsize; } } } } static void fwd_sums(double *f) { int stepsize,stage,curptr,nthreads,thread,step,nsteps; for(stage=dct_NxN_log2N-1; stage >=1; stage--){ nthreads = 1<<(stage-1); stepsize = nthreads<<1; nsteps = (1<<(dct_NxN_log2N-stage)) - 1; for(thread=1; thread<=nthreads; thread++){ curptr=nthreads +thread-1; for(step=1; step<=nsteps; step++){ f[curptr] += f[curptr+stepsize]; curptr += stepsize; } } } } // 打乱 static void scramble(double *f,int len){ int i,ii1,ii2; bitrev(f,len); bitrev(&f[0], len>>1); bitrev(&f[len>>1], len>>1); ii1=len-1; ii2=len>>1; for(i=0; i<(len>>2); i++){ SWAP(f[ii1], f[ii2]); ii1--; ii2++; } } // 恢复 static void unscramble(double *f,int len) { int i,ii1,ii2; ii1 = len-1; ii2 = len>>1; for(i=0; i<(len>>2); i++){ SWAP(f[ii1], f[ii2]); ii1--; ii2++; } bitrev(&f[0], len>>1); bitrev(&f[len>>1], len>>1); bitrev(f,len); } // 反转蝶形运算 static void inv_butterflies(double *f) { int stage,ii1,ii2,butterfly,ngroups,group,wingspan,increment,baseptr; double Cfac,T; for(stage=1; stage<=dct_NxN_log2N;stage++){ ngroups=1<<(dct_NxN_log2N-stage); wingspan=1<<(stage-1); increment=wingspan<<1; for(butterfly=1; butterfly<=wingspan; butterfly++){ Cfac = dct_NxN_costable[wingspan+butterfly-1]; baseptr=0; for(group=1; group<=ngroups; group++){ ii1=baseptr+butterfly-1; ii2=ii1+wingspan; T=Cfac * f[ii2]; f[ii2]=f[ii1]-T; f[ii1]=f[ii1]+T; baseptr += increment; } } } } // 前驱蝶形运算 static void fwd_butterflies(double *f) { int stage,ii1,ii2,butterfly,ngroups,group,wingspan,increment,baseptr; double Cfac,T; for(stage=dct_NxN_log2N; stage>=1;stage--){ ngroups=1<<(dct_NxN_log2N-stage); wingspan=1<<(stage-1); increment=wingspan<<1; for(butterfly=1; butterfly<=wingspan; butterfly++){ Cfac = dct_NxN_costable[wingspan+butterfly-1]; baseptr=0; for(group=1; group<=ngroups; group++){ ii1=baseptr+butterfly-1; ii2=ii1+wingspan; T= f[ii2]; f[ii2]=Cfac *(f[ii1]-T); f[ii1]=f[ii1]+T; baseptr += increment; } } } } // 逆尺度变换 static void ifct_noscale(double *f) { f[0] *= INVROOT2; inv_sums(f); bitrev(f,N); inv_butterflies(f); unscramble(f,N); } // 尺度变换 static void fct_noscale(double *f) { scramble(f,N); fwd_butterflies(f); bitrev(f,N); fwd_sums(f); f[0] *= INVROOT2; } // N×N的DCT void fdct_NxN(gray **pixels, double **dcts) { int u,v; double two_over_sqrtncolsnrows = 2.0/sqrt((double) N*M); for (u=0; u < N; u++) for (v=0; v < M; v++) dcts[u][v] = ((int) pixels[u][v] - 128); for (u=0; u<=M-1; u++){ fct_noscale(dcts[u]); } for (v=0; v<=N-1; v++){ for (u=0; u<=M-1; u++){ dct_NxN_tmp[u] = dcts[u][v]; } fct_noscale(dct_NxN_tmp); for (u=0; u<=M-1; u++){ dcts[u][v] = dct_NxN_tmp[u]*two_over_sqrtncolsnrows; } } } // N×N的逆DCT void idct_NxN(double **dcts, gray **pixels) { int u,v; double two_over_sqrtncolsnrows = 2.0/sqrt((double) N*M); double **tmp; tmp = alloc_coeffs(N, N); for (u=0;u<N;u++) for (v=0;v<M;v++) tmp[u][v] = dcts[u][v]; for (u=0; u<=M-1; u++){ ifct_noscale(tmp[u]); } for (v=0; v<=N-1; v++){ for (u=0; u<=M-1; u++){ dct_NxN_tmp[u] = tmp[u][v]; } ifct_noscale(dct_NxN_tmp); for (u=0; u<=M-1; u++){ tmp[u][v] = dct_NxN_tmp[u]*two_over_sqrtncolsnrows; } } for (u=0;u<N;u++) for (v=0;v<M;v++) pixels[u][v] = PIXELRANGE(tmp[u][v] + 128.5); free(tmp); } // N×N的inplace DCT void fdct_inplace_NxN(double **coeffs) { int u,v; double two_over_sqrtncolsnrows = 2.0/sqrt((double) N*M); for (u=0; u<=M-1; u++) fct_noscale(coeffs[u]); for (v=0; v<=N-1; v++){ for (u=0; u<=M-1; u++) dct_NxN_tmp[u] = coeffs[u][v]; fct_noscale(dct_NxN_tmp); for (u=0; u<=M-1; u++) coeffs[u][v] = dct_NxN_tmp[u]*two_over_sqrtncolsnrows; } } // N×N的inplace 逆DCT void idct_inplace_NxN(double **coeffs) { int u,v; double two_over_sqrtncolsnrows = 2.0/sqrt((double) N*M); for (u=0; u<=M-1; u++) ifct_noscale(coeffs[u]); for (v=0; v<=N-1; v++) { for (u=0; u<=M-1; u++) dct_NxN_tmp[u] = coeffs[u][v]; ifct_noscale(dct_NxN_tmp); for (u=0; u<=M-1; u++) coeffs[u][v] = dct_NxN_tmp[u]*two_over_sqrtncolsnrows; } } double **dct_NxM_co