void* iswt_2d(vector &swtop,int J, string nm, vector > &sig,int row, int col) { // 2D Inverse Stationary Wavelet Transform // Works In Conjunction with C++ Wavelet Library http://code.google.com/p/wavelet1d/ // (C) Rafat Hussain under GNU GPL v3 // See Project page http://code.google.com/p/wavelet1d/ // swtop [1D Output vector from swt_2d] // J [Number of Levels] // nm [Wavelet Name] // sig [Output] // row & col [Rows and cols of sig] // I have not fully tested this code. // References : // 1. Nason, G.P. and Silverman, B.W. (1995) The stationary wavelet transform and some statistical applications. // In Antoniadis, A. and Oppenheim, G. (eds) Lecture Notes in Statistics, 103, 281--300. // 2. Pesquet, J.C.; H. Krim, H. Carfatan (1996), "Time-invariant orthonormal wavelet representations," // http://pdf.aminer.org/000/318/848/frequency_shift_invariant_orthonormal_wavelet_packet_representations.pdf int rows_n =row; int cols_n = col; vector lp1,hp1,lp2,hp2; filtcoef(nm,lp1,hp1,lp2,hp2); vector low_pass = lp2; vector high_pass = hp2; int lf = low_pass.size(); int sum =0; vector > cLL(rows_n, vector(cols_n)); for (int iter=J; iter > 0; iter--) { vector > cLH(rows_n, vector(cols_n)); vector > cHL(rows_n, vector(cols_n)); vector > cHH(rows_n, vector(cols_n)); if (iter == J) { for (int i = 0 ; i < rows_n; i++ ){ for (int j = 0; j < cols_n; j++){ cLL[i][j] = swtop[sum+ i * cols_n + j]; cHL[i][j] = swtop[sum+ rows_n * cols_n+ i * cols_n + j]; cLH[i][j] = swtop[sum+ 2 * rows_n * cols_n + i * cols_n + j]; cHH[i][j] = swtop[sum+ 3* rows_n * cols_n + i * cols_n + j]; } } sum+= 4 * rows_n * cols_n; } else { for (int i = 0 ; i < rows_n; i++ ){ for (int j = 0; j < cols_n; j++){ cHL[i][j] = swtop[sum+ i * cols_n + j]; cLH[i][j] = swtop[sum+ rows_n * cols_n + i * cols_n + j]; cHH[i][j] = swtop[sum+ 2* rows_n * cols_n + i * cols_n + j]; } } sum+= 3 * rows_n * cols_n; } int value = (int) pow(2.0,(double) (iter-1)); for (int count = 0; count < value; count++) { vector > cLL1(rows_n/value, vector(cols_n/value)); vector > cLH1(rows_n/value, vector(cols_n/value)); vector > cHL1(rows_n/value, vector(cols_n/value)); vector > cHH1(rows_n/value, vector(cols_n/value)); int row1 = 0; int col1 = 0; for (int index_r = count; index_r < rows_n; index_r+=value){ for (int index_c = count; index_c < cols_n; index_c+=value){ cLL1[row1][col1]=cLL[index_r][index_c]; cLH1[row1][col1]=cLH[index_r][index_c]; cHL1[row1][col1]=cHL[index_r][index_c]; cHH1[row1][col1]=cHH[index_r][index_c]; col1++; } col1 = 0; row1++; } //Implementing shift =0 unsigned int len_c = cLL1[0].size(); unsigned int len_r = cLL1.size(); vector > cL(rows_n,vector(cols_n)); vector > cH(rows_n,vector(cols_n)); vector > cLL2(len_r/2, vector(len_c/2)); vector > cLH2(len_r/2, vector(len_c/2)); vector > cHL2(len_r/2, vector(len_c/2)); vector > cHH2(len_r/2, vector(len_c/2)); row1=0; col1=0; for (int index_r = 0; index_r < len_r; index_r+=2){ for (unsigned int index_c = 0; index_c < len_c; index_c+=2){ cLL2[row1][col1]=cLL1[index_r][index_c]; cLH2[row1][col1]=cLH1[index_r][index_c]; cHL2[row1][col1]=cHL1[index_r][index_c]; cHH2[row1][col1]=cHH1[index_r][index_c]; col1++; } col1 = 0; row1++; } vector > cLLu(len_r, vector(len_c)); vector > cLHu(len_r, vector(len_c)); vector > cHLu(len_r, vector(len_c)); vector > cHHu(len_r, vector(len_c)); int row_u = cLL2.size(); int col_u = cLL2[0].size(); upsamp2(cLL2,cLLu,2,2); upsamp2(cLH2,cLHu,2,2); upsamp2(cHL2,cHLu,2,2); upsamp2(cHH2,cHHu,2,2); //cLL2.clear();cLH2.clear();cHL2.clear();cHH2.clear(); vector > cLL00(len_r+lf, vector(len_c + lf)); vector > cLH00(len_r+lf, vector(len_c + lf)); vector > cHL00(len_r+lf, vector(len_c + lf)); vector > cHH00(len_r+lf, vector(len_c + lf)); per_ext2d(cLLu,cLL00,lf/2); per_ext2d(cLHu,cLH00,lf/2); per_ext2d(cHLu,cHL00,lf/2); per_ext2d(cHHu,cHH00,lf/2); //cLLu.clear();cLHu.clear();cHLu.clear();cHHu.clear(); //Shift By 1 int U = 2; // Upsampling Factor vector > cLL3(len_r/2, vector(len_c/2)); vector > cLH3(len_r/2, vector(len_c/2)); vector > cHL3(len_r/2, vector(len_c/2)); vector > cHH3(len_r/2, vector(len_c/2)); col1 = 0; row1 = 0; for (int index_r = 1; index_r < len_r; index_r+=2){ for (unsigned int index_c = 1; index_c < len_c; index_c+=2){ cLL3[row1][col1]=cLL1[index_r][index_c]; cLH3[row1][col1]=cLH1[index_r][index_c]; cHL3[row1][col1]=cHL1[index_r][index_c]; cHH3[row1][col1]=cHH1[index_r][index_c]; col1++; } col1 = 0; row1++; } vector > cLLu2(len_r, vector(len_c)); vector > cLHu2(len_r, vector(len_c)); vector > cHLu2(len_r, vector(len_c)); vector > cHHu2(len_r, vector(len_c)); row_u = cLL3.size(); col_u = cLL3[0].size(); upsamp2(cLL3,cLLu2,2,2); upsamp2(cLH3,cLHu2,2,2); upsamp2(cHL3,cHLu2,2,2); upsamp2(cHH3,cHHu2,2,2); //cLL3.clear();cLH3.clear();cHL3.clear();cHH3.clear(); cout << "STAGE 6" << endl; vector > cLL01(len_r+lf, vector(len_c + lf)); vector > cLH01(len_r+lf, vector(len_c + lf)); vector > cHL01(len_r+lf, vector(len_c + lf)); vector > cHH01(len_r+lf, vector(len_c + lf)); per_ext2d(cLLu2,cLL01,lf/2); per_ext2d(cLHu2,cLH01,lf/2); per_ext2d(cHLu2,cHL01,lf/2); per_ext2d(cHHu2,cHH01,lf/2); cLLu2.clear();cLHu2.clear();cHLu2.clear();cHHu2.clear(); int rowsLL = cLL01.size(); int colsLL = cLL01[0].size(); vector > cLLt00(rowsLL + lf-1, vector(colsLL)); vector > cLHt00(rowsLL + lf-1, vector(colsLL)); vector > cHLt00(rowsLL + lf-1, vector(colsLL)); vector > cHHt00(rowsLL + lf-1, vector(colsLL)); vector > cLLt10(rowsLL + lf-1, vector(colsLL)); vector > cLHt10(rowsLL + lf-1, vector(colsLL)); vector > cHLt10(rowsLL + lf-1, vector(colsLL)); vector > cHHt10(rowsLL + lf-1, vector(colsLL)); vector > cLLt01(rowsLL + lf-1, vector(colsLL + lf-1)); vector > cLHt01(rowsLL + lf-1, vector(colsLL + lf-1)); vector > cHLt01(rowsLL + lf-1, vector(colsLL + lf-1)); vector > cHHt01(rowsLL + lf-1, vector(colsLL + lf-1)); vector > cLLt11(rowsLL + lf-1, vector(colsLL + lf-1)); vector > cLHt11(rowsLL + lf-1, vector(colsLL + lf-1)); vector > cHLt11(rowsLL + lf-1, vector(colsLL + lf-1)); vector > cHHt11(rowsLL + lf-1, vector(colsLL + lf-1)); vector > cLLt0(rowsLL + lf-1, vector(colsLL + lf-1)); vector > cLLt1(rowsLL + lf-1, vector(colsLL + lf-1)); vector > oupt0(len_r, vector(len_c)); vector > oupt1(len_r, vector(len_c)); vector > oupt(len_r, vector(len_c)); for (int j=0; j< colsLL;j++) { vector sigLL0; vector sigLH0; vector sigHL0; vector sigHH0; for (int i=0; i < rowsLL; i++) { double temp01 = cLL00[i][j]; double temp02 = cLH00[i][j]; double temp03 = cHL00[i][j]; double temp04 = cHH00[i][j]; sigLL0.push_back(temp01); sigLH0.push_back(temp02); sigHL0.push_back(temp03); sigHH0.push_back(temp04); } // First Convolution Step for LL,LH,HL,HH vector X0LL; convfft(sigLL0, low_pass, X0LL); vector X0LH; convfft(sigLH0, low_pass, X0LH); vector X0HL; convfft(sigHL0, high_pass, X0HL); vector X0HH; convfft(sigHH0, high_pass, X0HH); int lent=X0LL.size(); for (int i=0; i < lent; i++) { cLLt00[i][j]=X0LL[i]; cLHt00[i][j]=X0LH[i]; cHLt00[i][j]=X0HL[i]; cHHt00[i][j]=X0HH[i]; } } for (int i=0; i< rowsLL+lf-1;i++) { vector sigLL0; vector sigLH0; vector sigHL0; vector sigHH0; for (int j=0; j < colsLL; j++) { double temp01 = cLLt00[i][j]; double temp02 = cLHt00[i][j]; double temp03 = cHLt00[i][j]; double temp04 = cHHt00[i][j]; sigLL0.push_back(temp01); sigLH0.push_back(temp02); sigHL0.push_back(temp03); sigHH0.push_back(temp04); } // Second Convolution Step for LL,LH,HL,HH vector X0LL; convfft(sigLL0, low_pass, X0LL); vector X0LH; convfft(sigLH0, high_pass, X0LH); vector X0HL; convfft(sigHL0, low_pass, X0HL); vector X0HH; convfft(sigHH0, high_pass, X0HH); int lent=X0LL.size(); for (int j=0; j < lent; j++) { cLLt01[i][j]=X0LL[j]; cLHt01[i][j]=X0LH[j]; cHLt01[i][j]=X0HL[j]; cHHt01[i][j]=X0HH[j]; cLLt0[i][j]=X0LL[j]+X0LH[j]+X0HL[j]+X0HH[j]; } } for (int j=0; j< colsLL;j++) { vector sigLL0; vector sigLH0; vector sigHL0; vector sigHH0; for (int i=0; i < rowsLL; i++) { double temp01 = cLL01[i][j]; double temp02 = cLH01[i][j]; double temp03 = cHL01[i][j]; double temp04 = cHH01[i][j]; sigLL0.push_back(temp01); sigLH0.push_back(temp02); sigHL0.push_back(temp03); sigHH0.push_back(temp04); } // First Convolution Step for LL,LH,HL,HH vector X0LL; convfft(sigLL0, low_pass, X0LL); vector X0LH; convfft(sigLH0, low_pass, X0LH); vector X0HL; convfft(sigHL0, high_pass, X0HL); vector X0HH; convfft(sigHH0, high_pass, X0HH); int lent=X0LL.size(); for (int i=0; i < lent; i++) { cLLt10[i][j]=X0LL[i]; cLHt10[i][j]=X0LH[i]; cHLt10[i][j]=X0HL[i]; cHHt10[i][j]=X0HH[i]; } } for (int i=0; i< rowsLL + lf-1;i++) { vector sigLL0; vector sigLH0; vector sigHL0; vector sigHH0; for (int j=0; j < colsLL; j++) { double temp01 = cLLt10[i][j]; double temp02 = cLHt10[i][j]; double temp03 = cHLt10[i][j]; double temp04 = cHHt10[i][j]; sigLL0.push_back(temp01); sigLH0.push_back(temp02); sigHL0.push_back(temp03); sigHH0.push_back(temp04); } // Second Convolution Step for LL,LH,HL,HH vector X0LL; convfft(sigLL0, low_pass, X0LL); vector X0LH; convfft(sigLH0, high_pass, X0LH); vector X0HL; convfft(sigHL0, low_pass, X0HL); vector X0HH; convfft(sigHH0, high_pass, X0HH); int lent=X0LL.size(); for (int j=0; j < lent; j++) { cLLt11[i][j]=X0LL[j]; cLHt11[i][j]=X0LH[j]; cHLt11[i][j]=X0HL[j]; cHHt11[i][j]=X0HH[j]; cLLt1[i][j]=X0LL[j]+X0LH[j]+X0HL[j]+X0HH[j]; } } for (int i=0; i < len_r;i++) { for (int j=0; j < len_c;j++) { oupt0[i][j]=cLLt0[lf+i-1][lf+j-1]; oupt1[i][j]=cLLt1[lf+i-1][lf+j-1]; } } circshift2d(oupt1,-1,-1); for (int i=0; i < len_r;i++) { for (int j=0; j < len_c;j++) { double temp=oupt0[i][j]+oupt1[i][j]; temp=temp/2.0; oupt[i][j]=temp; } } row1 = 0; col1 = 0; cout << rows_n << cols_n << " " << oupt.size() << oupt[0].size() << endl; for (int index_r = count; index_r < rows_n; index_r+=value){ for (int index_c = count; index_c < cols_n; index_c+=value){ cLL[index_r][index_c]=oupt[row1][col1]; col1++; } col1 = 0; row1++; } cout << "STAGE 14" << endl; } if (iter == 1) { sig=cLL; } } return 0; }