/* nag_dwt_3d (c09fac) Example Program. * * Copyright 2013 Numerical Algorithms Group. * * Mark 24, 2013. */ #include #include #include #include #include #define A(I,J,K) a[I-1 + (J-1)* lda + (K-1)* lda * sda] #define B(I,J,K) b[I-1 + (J-1)* ldb + (K-1)* ldb * sdb] #define D(I,J,K) d[I-1 + (J-1)* nwcm + (K-1)* nwcm * nwcn] int main(void) { /* Scalars */ Integer exit_status = 0, zero = 0; Integer cindex, i, j, k, lda, ldb, lenc; Integer m, n, fr, nf, nwcfr, nwcm, nwcn, nwct, nwl, sda, sdb; /* Arrays */ char mode[25], wavnam[25]; double *a = 0, *b = 0, *c = 0, *d = 0; Integer icomm[260]; /* Nag Types */ Nag_Wavelet wavnamenum; Nag_WaveletMode modenum; Nag_MatrixType matrix = Nag_GeneralMatrix; Nag_OrderType order = Nag_ColMajor; Nag_DiagType diag = Nag_NonUnitDiag; NagError fail; INIT_FAIL(fail); printf("nag_dwt_3d (c09fac) Example Program Results\n\n"); /* Skip heading in data file and read problem parameters */ scanf("%*[^\n] %ld%ld%ld%*[^\n]", &m, &n, &fr); lda = m; ldb = m; sda = n; sdb = n; scanf("%25s%25s%*[^\n]\n", wavnam, mode); if (!(a = NAG_ALLOC((lda)*(sda)*(fr), double)) || !(b = NAG_ALLOC((ldb)*(sdb)*(fr), double))) { printf("Allocation failure\n"); exit_status = 1; goto END; } printf("Parameters read from file :: \n"); printf("DWT :: Wavelet : %s\n", wavnam); printf(" End mode : %s\n", mode); printf(" m : %4ld\n", m); printf(" n : %4ld\n", n); printf(" fr : %4ld\n\n", fr); /* nag_enum_name_to_value (x04nac). * Converts NAG enum member name to value */ wavnamenum = (Nag_Wavelet) nag_enum_name_to_value(wavnam); modenum = (Nag_WaveletMode) nag_enum_name_to_value(mode); /* Read data array*/ for (k=1; k<=fr; k++) { for (i=1; i<=m; i++) { for (j=1; j<=n; j++) scanf("%lf", &A(i, j, k)); } scanf("%*[^\n] "); } /* Print out the input data */ printf("Input Data :\n"); for (k=1; k<=fr; k++) { /* nag_gen_real_mat_print_comp (x04cbc). * Prints out a matrix. */ nag_gen_real_mat_print_comp(order, matrix, diag, m, n, &A(1, 1, k), lda, "%8.4f"," ", Nag_NoLabels, 0, Nag_NoLabels, 0, 80, 0, 0, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_gen_real_mat_print_comp (x04cbc).\n%s\n", fail.message); exit_status = 2; goto END; } printf("\n"); } /* nag_wfilt_3d (c09acc). * Three-dimensional wavelet filter initialization */ nag_wfilt_3d(wavnamenum, Nag_SingleLevel, modenum, m, n, fr, &nwl, &nf, &nwct, &nwcn, &nwcfr, icomm, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_wfilt_3d (c09acc).\n%s\n", fail.message); exit_status = 3; goto END; } /* Calculate the number of wavelet coefficients in * the first dimension, nwcm. */ nwcm = nwct/(8*nwcn*nwcfr); lenc = nwct; /* Allocate space for the coefficients array, C */ if (!(c = NAG_ALLOC((lenc), double))) { printf("Allocation failure\n"); exit_status = 4; goto END; } /* nag_dwt_3d (c09fac). * Three-dimensional discrete wavelet transform */ nag_dwt_3d(m, n, fr, a, lda, sda, lenc, c, icomm, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_dwt_3d (c09fac).\n%s\n", fail.message); exit_status = 5; goto END; } /* Allocate space for extraction of coefficients of a single type */ if (!(d = NAG_ALLOC((nwcm)*(nwcn)*(nwcfr), double))) { printf("Allocation failure\n"); exit_status = 6; goto END; } for (cindex=0; cindex<=7; cindex++) { /* Use the extraction routine c09fyc to retrieve the required * coefficients. */ /* nag_wav_3d_coeff_ext (c09fyc). * Extract the nominated coefficients. */ nag_wav_3d_coeff_ext(zero,cindex,lenc,c,d,nwcm,nwcn,icomm,&fail); if (fail.code != NE_NOERROR) { printf("Error from nag_wav_3d_coeff_ext (c09fyc).\n%s\n", fail.message); exit_status = 7; goto END; } /* Print out the extracted coefficients */ switch (cindex) { case 0: printf("Approximation coefficients (LLL)\n"); break; case 1: printf("Detail coefficients (LLH)\n"); break; case 2: printf("Detail coefficients (LHL)\n"); break; case 3: printf("Detail coefficients (LHH)\n"); break; case 4: printf("Detail coefficients (HLL)\n"); break; case 5: printf("Detail coefficients (HLH)\n"); break; case 6: printf("Detail coefficients (HHL)\n"); break; case 7: printf("Detail coefficients (HHH)\n"); break; } for (i=1; i<=nwcm; i++) { if (i==1) { printf("Coefficients "); for (k=1; k<=nwcfr; k++) { printf("Frame %4"NAG_IFMT, k); for (j=1; j<=9*nwcn-8; j++) printf(" "); } printf("\n"); } for (k=1; k<=nwcfr; k++) { if (k==1 && i==1) printf("%5ld%8s", cindex, " "); else if (k==1) printf("%13s"," "); else printf("%2s"," "); for (j=1; j<=nwcn; j++) { printf("%8.4f ",D(i,j,k)); } } printf("\n"); } printf("\n"); } /* nag_idwt_3d (c09fbc). * Three-dimensional inverse discrete wavelet transform */ nag_idwt_3d(m, n, fr, lenc, c, b, ldb, sdb, icomm, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_idwt_3d (c09fbc).\n%s\n", fail.message); exit_status = 8; goto END; } printf("Output Data :\n"); for (k=1; k<=fr; k++) { /* nag_gen_real_mat_print_comp (x04cbc). * Prints out a matrix. */ nag_gen_real_mat_print_comp(order, matrix, diag, m, n, &B(1, 1, k), ldb, "%8.4f"," ", Nag_NoLabels, 0, Nag_NoLabels, 0, 80, 0, 0, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_gen_real_mat_print_comp (x04cbc).\n%s\n", fail.message); exit_status = 9; goto END; } printf("\n"); } END: NAG_FREE(a); NAG_FREE(b); NAG_FREE(c); NAG_FREE(d); return exit_status; }