1 /*
2 * -------------------------------------------------------------------------
3 * dwt2d.c -- 2-D signal decomposition and reconstruction
4 * SWT - Scilab wavelet toolbox
5 * Copyright (C) 2005-2006 Roger Liu
6 * Copyright (C) 20010-2012 Holger Nahrstaedt
7 *
8 * This program is free software; you can redistribute it and/or modify
9 * it under the terms of the GNU General Public License as published by
10 * the Free Software Foundation; either version 2 of the License, or
11 * (at your option) any later version.
12 *
13 * This program is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 * GNU General Public License for more details.
17 *
18 * You should have received a copy of the GNU General Public License
19 * along with this program; if not, write to the Free Software
20 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
21 * -------------------------------------------------------------------------
22 */
23
24 #include "swtlib.h"
25
26 void
dwt2D(double * matrixIn,int matrixInRow,int matrixInCol,double * matrixOutApprox,double * matrixOutColDetail,double * matrixOutRowDetail,double * matrixOutDetail,int matrixOutRow,int matrixOutCol,double * lowDe,double * hiDe,int filterLen,extend_method extMethod)27 dwt2D (double *matrixIn, int matrixInRow, int matrixInCol,
28 double *matrixOutApprox, double *matrixOutColDetail,
29 double *matrixOutRowDetail, double *matrixOutDetail,
30 int matrixOutRow, int matrixOutCol, double *lowDe,
31 double *hiDe, int filterLen, extend_method extMethod)
32 {
33 int row, col;
34 int filterOutLength;
35 int matrixOutRowM, matrixOutColM, matrixInRowM, matrixInColM;
36 //double temp;
37 double *matrixInPre, *matrixInRo;
38 double *pMatrixOutApproxPre, *pMatrixOutDetailPre;
39 double *matrixOutApproxPre, *matrixOutDetailPre;
40 double *matrixOutApproxM, *matrixOutColDetailM;
41 double *matrixOutRowDetailM, *matrixOutDetailM;
42 char c='b';
43
44 /* extension */
45 //matrixInRowM = matrixInRow + 2 * (filterLen - 1);
46 //matrixInColM = matrixInCol + 2 * (filterLen - 1);
47 matrixInRowM = matrixInRow + 2 * filterLen;
48 matrixInColM = matrixInCol + 2 * filterLen;
49 if ((extMethod==PER)&&(matrixInRow%2!=0))
50 matrixInRowM = matrixInRow + 2 * filterLen + 1;
51 if ((extMethod==PER)&&(matrixInCol%2!=0))
52 matrixInColM = matrixInCol + 2 * filterLen + 1;
53
54 matrixInRo = malloc (matrixInRowM * matrixInColM * sizeof (double));
55 matrixInPre = malloc(matrixInRowM * matrixInColM * sizeof (double));
56 wextend_2D (matrixIn, matrixInRow, matrixInCol, matrixInRo,
57 matrixInRowM, matrixInColM, extMethod, &c, &c);
58 matrix_tran (matrixInRo, matrixInColM, matrixInRowM,
59 matrixInPre, matrixInRowM, matrixInColM);
60 free (matrixInRo);
61
62 /* Row DWT */
63
64 filterOutLength = matrixInColM + filterLen - 1;
65
66 matrixOutColM = filterOutLength/2;
67 filterOutLength = matrixInRowM + filterLen - 1;
68 //if ((extMethod==PER)&&(matrixInRow%2!=0))
69 // filterOutLength =
70 matrixOutRowM = filterOutLength/2;
71
72 matrixOutApproxPre = malloc (matrixInRowM *
73 matrixOutColM * sizeof (double));
74 matrixOutDetailPre = malloc (matrixInRowM *
75 matrixOutColM * sizeof (double));
76 filterOutLength = matrixInColM + filterLen - 1;
77 for (row = 0; row < matrixInRowM; row++)
78 {
79 dwt_no_extension ((matrixInPre + row * matrixInColM), matrixInColM, lowDe,
80 hiDe, filterLen,
81 (matrixOutApproxPre + row * matrixOutColM),
82 (matrixOutDetailPre + row * matrixOutColM),
83 matrixOutColM);//, extMethod);
84 }
85
86 free (matrixInPre);
87
88 /* transpose */
89 pMatrixOutApproxPre = malloc (matrixInRowM *
90 matrixOutColM * sizeof (double));
91 matrix_tran (matrixOutApproxPre, matrixInRowM, matrixOutColM,
92 pMatrixOutApproxPre, matrixInRowM, matrixOutColM);
93 free (matrixOutApproxPre);
94
95 pMatrixOutDetailPre = malloc (matrixInRowM *
96 matrixOutColM * sizeof (double));
97 matrix_tran (matrixOutDetailPre, matrixInRowM, matrixOutColM,
98 pMatrixOutDetailPre, matrixInRowM, matrixOutColM);
99
100 free (matrixOutDetailPre);
101
102 /* Col DWT */
103 filterOutLength = matrixInRowM + filterLen - 1;
104 matrixOutApproxM = malloc (matrixOutRowM * matrixOutColM
105 * sizeof (double));
106 matrixOutColDetailM = malloc (matrixOutRowM * matrixOutColM
107 * sizeof (double));
108
109 for (col = 0; col < matrixOutColM; col++)
110 {
111 dwt_no_extension ((pMatrixOutApproxPre + col * matrixInRowM), matrixInRowM,
112 lowDe, hiDe, filterLen,
113 (matrixOutApproxM + col * matrixOutRowM),
114 (matrixOutColDetailM + col * matrixOutRowM),
115 matrixOutRowM);//, extMethod);
116 }
117 free (pMatrixOutApproxPre);
118
119 wkeep_2D_center (matrixOutApproxM, matrixOutRowM, matrixOutColM,
120 matrixOutApprox, matrixOutRow, matrixOutCol);
121 free (matrixOutApproxM);
122 wkeep_2D_center (matrixOutColDetailM, matrixOutRowM, matrixOutColM,
123 matrixOutColDetail, matrixOutRow, matrixOutCol);
124 free (matrixOutColDetailM);
125
126 filterOutLength = matrixInRowM + filterLen - 1;
127 matrixOutRowDetailM = malloc (matrixOutRowM * matrixOutColM
128 * sizeof (double));
129 matrixOutDetailM = malloc (matrixOutRowM * matrixOutColM * sizeof (double));
130 for (col = 0; col < matrixOutColM; col++)
131 {
132 dwt_no_extension ((pMatrixOutDetailPre + col * matrixInRowM), matrixInRowM,
133 lowDe, hiDe, filterLen,
134 (matrixOutRowDetailM + col * matrixOutRowM),
135 (matrixOutDetailM + col * matrixOutRowM),
136 matrixOutRowM);//, extMethod);
137 }
138 free (pMatrixOutDetailPre);
139 wkeep_2D_center (matrixOutRowDetailM, matrixOutRowM, matrixOutColM,
140 matrixOutRowDetail, matrixOutRow, matrixOutCol);
141 free (matrixOutRowDetailM);
142 wkeep_2D_center (matrixOutDetailM, matrixOutRowM, matrixOutColM,
143 matrixOutDetail, matrixOutRow, matrixOutCol);
144 free (matrixOutDetailM);
145
146 return;
147 }
148
149
150 void
dwt2D_neo(double * matrixIn,int matrixInRow,int matrixInCol,double * matrixOutApprox,double * matrixOutColDetail,double * matrixOutRowDetail,double * matrixOutDetail,int matrixOutRow,int matrixOutCol,double * lowDe,double * hiDe,int filterLen,extend_method extMethod)151 dwt2D_neo (double *matrixIn, int matrixInRow, int matrixInCol,
152 double *matrixOutApprox, double *matrixOutColDetail,
153 double *matrixOutRowDetail, double *matrixOutDetail,
154 int matrixOutRow, int matrixOutCol, double *lowDe,
155 double *hiDe, int filterLen, extend_method extMethod)
156 {
157 int row, col;
158 int filterOutLength;
159 int matrixOutRowM, matrixOutColM, matrixInRowM, matrixInColM;
160 int matrixOutRowT, matrixOutColT;
161 double *matrixInPre, *matrixInRo;
162 double *pMatrixOutApproxPre, *pMatrixOutDetailPre;
163 double *matrixOutApproxPre, *matrixOutDetailPre;
164 double *matrixOutApproxM, *matrixOutColDetailM;
165 double *matrixOutRowDetailM, *matrixOutDetailM;
166 double *matrixOutApproxT, *matrixOutColDetailT;
167 double *matrixOutRowDetailT, *matrixOutDetailT;
168 char c='b';
169
170 /* extension */
171 matrixInRowM = matrixInRow + 2 * filterLen;
172 matrixInColM = matrixInCol + 2 * filterLen;
173 if ((extMethod==PER)&&(matrixInRow%2!=0))
174 matrixInRowM = matrixInRow + 2 * filterLen + 1;
175 if ((extMethod==PER)&&(matrixInCol%2!=0))
176 matrixInColM = matrixInCol + 2 * filterLen + 1;
177
178 matrixInRo = malloc (matrixInRowM * matrixInColM * sizeof (double));
179 matrixInPre = malloc(matrixInRowM * matrixInColM * sizeof (double));
180 wextend_2D (matrixIn, matrixInRow, matrixInCol, matrixInRo,
181 matrixInRowM, matrixInColM, extMethod, &c, &c);
182 matrix_tran (matrixInRo, matrixInColM, matrixInRowM,
183 matrixInPre, matrixInRowM, matrixInColM);
184 free (matrixInRo);
185
186 /* Row DWT */
187
188 filterOutLength = matrixInColM + filterLen - 1;
189 matrixOutColM = filterOutLength;
190 filterOutLength = matrixInRowM + filterLen - 1;
191 matrixOutRowM = filterOutLength;
192
193 matrixOutApproxPre = malloc (matrixInRowM *
194 matrixOutColM * sizeof (double));
195 matrixOutDetailPre = malloc (matrixInRowM *
196 matrixOutColM * sizeof (double));
197 filterOutLength = matrixInColM + filterLen - 1;
198 for (row = 0; row < matrixInRowM; row++)
199 {
200 dwt_conv ((matrixInPre + row * matrixInColM), matrixInColM, lowDe,
201 hiDe, filterLen,
202 (matrixOutApproxPre + row * matrixOutColM),
203 (matrixOutDetailPre + row * matrixOutColM),
204 matrixOutColM);
205 }
206
207 free (matrixInPre);
208
209 /* transpose */
210 pMatrixOutApproxPre = malloc (matrixInRowM *
211 matrixOutColM * sizeof (double));
212 matrix_tran (matrixOutApproxPre, matrixInRowM, matrixOutColM,
213 pMatrixOutApproxPre, matrixInRowM, matrixOutColM);
214 free (matrixOutApproxPre);
215
216 pMatrixOutDetailPre = malloc (matrixInRowM *
217 matrixOutColM * sizeof (double));
218 matrix_tran (matrixOutDetailPre, matrixInRowM, matrixOutColM,
219 pMatrixOutDetailPre, matrixInRowM, matrixOutColM);
220
221 free (matrixOutDetailPre);
222
223 /* Col DWT */
224 filterOutLength = matrixInRowM + filterLen - 1;
225 matrixOutApproxM = malloc (matrixOutRowM * matrixOutColM
226 * sizeof (double));
227 matrixOutColDetailM = malloc (matrixOutRowM * matrixOutColM
228 * sizeof (double));
229
230 //printf("extmethod %d",extMethod);
231 for (col = 0; col < matrixOutColM; col++)
232 {
233 dwt_conv ((pMatrixOutApproxPre + col * matrixInRowM), matrixInRowM,
234 lowDe, hiDe, filterLen,
235 (matrixOutApproxM + col * matrixOutRowM),
236 (matrixOutColDetailM + col * matrixOutRowM),
237 matrixOutRowM);//, extMethod);
238 }
239 free (pMatrixOutApproxPre);
240
241 matrixOutRowT = matrixInRow + filterLen - 1;
242 matrixOutColT = matrixInCol + filterLen - 1;
243 if ((extMethod==PER)&&(matrixInRow%2!=0))
244 matrixOutRowT = matrixInRow + 1;
245 if ((extMethod==PER)&&(matrixInCol%2!=0))
246 matrixOutColT = matrixInCol + 1;
247 if ((extMethod==PER)&&(matrixInRow%2==0))
248 matrixOutRowT = matrixInRow;
249 if ((extMethod==PER)&&(matrixInCol%2==0))
250 matrixOutColT = matrixInCol;
251 matrixOutApproxT = malloc (matrixOutRowT*matrixOutColT*sizeof(double));
252 matrixOutColDetailT = malloc (matrixOutRowT*matrixOutColT*sizeof(double));
253
254 wkeep_2D_center (matrixOutApproxM, matrixOutRowM, matrixOutColM,
255 matrixOutApproxT, matrixOutRowT, matrixOutColT);
256 free (matrixOutApproxM);
257 dyaddown_2D_keep_even (matrixOutApproxT, matrixOutRowT,
258 matrixOutColT, matrixOutApprox,
259 matrixOutRow, matrixOutCol);
260 free(matrixOutApproxT);
261
262 wkeep_2D_center (matrixOutColDetailM, matrixOutRowM, matrixOutColM,
263 matrixOutColDetailT, matrixOutRowT, matrixOutColT);
264 free (matrixOutColDetailM);
265 dyaddown_2D_keep_even (matrixOutColDetailT, matrixOutRowT,
266 matrixOutColT, matrixOutColDetail,
267 matrixOutRow, matrixOutCol);
268 free(matrixOutColDetailT);
269
270
271
272
273 filterOutLength = matrixInRowM + filterLen - 1;
274 matrixOutRowDetailM = malloc (matrixOutRowM * matrixOutColM
275 * sizeof (double));
276 matrixOutDetailM = malloc (matrixOutRowM * matrixOutColM * sizeof (double));
277 for (col = 0; col < matrixOutColM; col++)
278 {
279 dwt_conv ((pMatrixOutDetailPre + col * matrixInRowM), matrixInRowM,
280 lowDe, hiDe, filterLen,
281 (matrixOutRowDetailM + col * matrixOutRowM),
282 (matrixOutDetailM + col * matrixOutRowM),
283 matrixOutRowM);
284 }
285 free (pMatrixOutDetailPre);
286 matrixOutRowDetailT = malloc (matrixOutRowT*matrixOutColT*sizeof(double));
287 matrixOutDetailT = malloc (matrixOutRowT*matrixOutColT*sizeof(double));
288
289
290 wkeep_2D_center (matrixOutRowDetailM, matrixOutRowM, matrixOutColM,
291 matrixOutRowDetailT, matrixOutRowT, matrixOutColT);
292 free (matrixOutRowDetailM);
293 dyaddown_2D_keep_even (matrixOutRowDetailT, matrixOutRowT,
294 matrixOutColT, matrixOutRowDetail,
295 matrixOutRow, matrixOutCol);
296 free(matrixOutRowDetailT);
297
298 wkeep_2D_center (matrixOutDetailM, matrixOutRowM, matrixOutColM,
299 matrixOutDetailT, matrixOutRowT, matrixOutColT);
300 free (matrixOutDetailM);
301 dyaddown_2D_keep_even (matrixOutDetailT, matrixOutRowT,
302 matrixOutColT, matrixOutDetail,
303 matrixOutRow, matrixOutCol);
304 free(matrixOutDetailT);
305
306 return;
307 }
308
309
310 void
idwt2D(double * matrixInApprox,double * matrixInColDetail,double * matrixInRowDetail,double * matrixInDetail,int matrixInRow,int matrixInCol,double * lowRe,double * hiRe,int filterLen,double * matrixOut,int matrixOutRow,int matrixOutCol,extend_method extMethod)311 idwt2D (double *matrixInApprox, double *matrixInColDetail,
312 double *matrixInRowDetail, double *matrixInDetail,
313 int matrixInRow, int matrixInCol, double *lowRe,
314 double *hiRe, int filterLen, double *matrixOut,
315 int matrixOutRow, int matrixOutCol, extend_method extMethod)
316 {
317 int row, col;
318 int filterOutLength;
319 int matrixInRowM, matrixInColM;
320 char c='b';
321 double *matrixOutApproxPre, *matrixOutDetailPre;
322 double *matrixOutApproxTemp, *matrixOutDetailTemp;
323 double *matrixOutPre;
324 double *matrixInApproxM, *matrixInColDetailM;
325 double *matrixInRowDetailM, *matrixInDetailM;
326
327 matrixInRowM = matrixInRow + 2 * (filterLen - 1);
328 matrixInColM = matrixInCol + 2 * (filterLen - 1);
329
330 /* Approx extension */
331 matrixInApproxM = malloc (matrixInRowM * matrixInColM * sizeof (double));
332 wextend_2D (matrixInApprox, matrixInRow, matrixInCol,
333 matrixInApproxM, matrixInRowM, matrixInColM,
334 extMethod, &c, &c);
335
336 matrixInColDetailM = malloc (matrixInRowM * matrixInColM * sizeof (double));
337 wextend_2D (matrixInColDetail, matrixInRow, matrixInCol,
338 matrixInColDetailM, matrixInRowM, matrixInColM,
339 extMethod, &c, &c);
340
341 matrixInRowDetailM = malloc (matrixInRowM * matrixInColM * sizeof (double));
342 wextend_2D (matrixInRowDetail, matrixInRow, matrixInCol,
343 matrixInRowDetailM, matrixInRowM, matrixInColM,
344 extMethod, &c, &c);
345
346 matrixInDetailM = malloc (matrixInRowM * matrixInColM * sizeof (double));
347
348 wextend_2D (matrixInDetail, matrixInRow, matrixInCol,
349 matrixInDetailM, matrixInRowM, matrixInColM,
350 extMethod, &c, &c);
351
352 filterOutLength = 2 * matrixInRowM + filterLen - 1;
353
354 /* Approx Calculation */
355 matrixOutApproxPre = malloc (matrixOutRow * matrixInColM * sizeof (double));
356 matrixOutApproxTemp =
357 malloc (matrixOutRow * matrixInColM * sizeof (double));
358 for (col = 0; col < matrixInColM; col++)
359 idwt_neo ((matrixInApproxM + col * matrixInRowM),
360 (matrixInColDetailM + col * matrixInRowM),
361 matrixInRowM, lowRe, hiRe, filterLen,
362 (matrixOutApproxPre + col * matrixOutRow),
363 matrixOutRow);
364 matrix_tran (matrixOutApproxPre, matrixInColM, matrixOutRow,
365 matrixOutApproxTemp, matrixInColM, matrixOutRow);
366 free (matrixOutApproxPre);
367 free (matrixInApproxM);
368 free (matrixInColDetailM);
369
370 /* Detail Calculation */
371 matrixOutDetailPre = malloc (matrixOutRow * matrixInColM * sizeof (double));
372 for (col = 0; col < matrixInColM; col++)
373 idwt_neo ((matrixInRowDetailM + col * matrixInRowM),
374 (matrixInDetailM + col * matrixInRowM),
375 matrixInRowM, lowRe, hiRe, filterLen,
376 (matrixOutDetailPre + col * matrixOutRow),
377 matrixOutRow);
378 matrixOutDetailTemp =
379 malloc (matrixOutRow * matrixInColM * sizeof (double));
380 matrix_tran (matrixOutDetailPre, matrixInColM, matrixOutRow,
381 matrixOutDetailTemp, matrixInColM, matrixOutRow);
382 free (matrixOutDetailPre);
383 free (matrixInRowDetailM);
384 free (matrixInDetailM);
385
386 /* Final Inverse Transform */
387 filterOutLength = 2 * matrixInColM + filterLen - 1;
388 matrixOutPre = malloc (matrixOutRow * matrixOutCol * sizeof (double));
389 for (row = 0; row < matrixOutRow; row++)
390 idwt_neo ((matrixOutApproxTemp + row * matrixInColM),
391 (matrixOutDetailTemp + row * matrixInColM),
392 matrixInColM, lowRe, hiRe, filterLen,
393 (matrixOutPre + row * matrixOutCol),
394 matrixOutCol);
395 free (matrixOutApproxTemp);
396 free (matrixOutDetailTemp);
397 matrix_tran (matrixOutPre, matrixOutRow, matrixOutCol, matrixOut,
398 matrixOutRow, matrixOutCol);
399 free (matrixOutPre);
400
401 return;
402 }
403
404
405 void
idwt2D_neo(double * matrixInApprox,double * matrixInColDetail,double * matrixInRowDetail,double * matrixInDetail,int matrixInRow,int matrixInCol,double * lowRe,double * hiRe,int filterLen,double * matrixOut,int matrixOutRow,int matrixOutCol)406 idwt2D_neo (double *matrixInApprox, double *matrixInColDetail,
407 double *matrixInRowDetail, double *matrixInDetail,
408 int matrixInRow, int matrixInCol, double *lowRe,
409 double *hiRe, int filterLen, double *matrixOut,
410 int matrixOutRow, int matrixOutCol)
411 {
412 int row, col;
413 int filterOutLength;
414 int matrixInRowM, matrixInColM;
415 char c='b';
416 double *matrixOutApproxPre, *matrixOutDetailPre;
417 double *matrixOutApproxTemp, *matrixOutDetailTemp;
418 double *matrixOutPre;
419 //double *matrixInApproxM, *matrixInColDetailM;
420 //double *matrixInRowDetailM, *matrixInDetailM;
421
422 matrixInRowM = matrixInRow;// + 2 * (filterLen - 1);
423 matrixInColM = matrixInCol;// + 2 * (filterLen - 1);
424
425 /* Approx extension */
426 //matrixInApproxM = malloc (matrixInRowM * matrixInColM * sizeof (double));
427 //wextend_2D (matrixInApprox, matrixInRow, matrixInCol,
428 // matrixInApproxM, matrixInRowM, matrixInColM,
429 // extMethod, &c, &c);
430
431 //matrixInColDetailM = malloc (matrixInRowM * matrixInColM * sizeof (double));
432 //wextend_2D (matrixInColDetail, matrixInRow, matrixInCol,
433 // matrixInColDetailM, matrixInRowM, matrixInColM,
434 // extMethod, &c, &c);
435
436 //matrixInRowDetailM = malloc (matrixInRowM * matrixInColM * sizeof (double));
437 //wextend_2D (matrixInRowDetail, matrixInRow, matrixInCol,
438 // matrixInRowDetailM, matrixInRowM, matrixInColM,
439 // extMethod, &c, &c);
440
441 //matrixInDetailM = malloc (matrixInRowM * matrixInColM * sizeof (double));
442
443 //wextend_2D (matrixInDetail, matrixInRow, matrixInCol,
444 // matrixInDetailM, matrixInRowM, matrixInColM,
445 // extMethod, &c, &c);
446
447 filterOutLength = 2 * matrixInRowM + filterLen - 1;
448
449 /* Approx Calculation */
450 matrixOutApproxPre = malloc (matrixOutRow * matrixInColM * sizeof (double));
451 matrixOutApproxTemp =
452 malloc (matrixOutRow * matrixInColM * sizeof (double));
453 for (col = 0; col < matrixInColM; col++)
454 idwt_neo ((matrixInApprox + col * matrixInRowM),
455 (matrixInColDetail + col * matrixInRowM),
456 matrixInRowM, lowRe, hiRe, filterLen,
457 (matrixOutApproxPre + col * matrixOutRow),
458 matrixOutRow);
459 matrix_tran (matrixOutApproxPre, matrixInColM, matrixOutRow,
460 matrixOutApproxTemp, matrixInColM, matrixOutRow);
461 free (matrixOutApproxPre);
462 //free (matrixInApproxM);
463 //free (matrixInColDetailM);
464
465 /* Detail Calculation */
466 matrixOutDetailPre = malloc (matrixOutRow * matrixInColM * sizeof (double));
467 for (col = 0; col < matrixInColM; col++)
468 idwt_neo ((matrixInRowDetail + col * matrixInRowM),
469 (matrixInDetail + col * matrixInRowM),
470 matrixInRowM, lowRe, hiRe, filterLen,
471 (matrixOutDetailPre + col * matrixOutRow),
472 matrixOutRow);
473 matrixOutDetailTemp =
474 malloc (matrixOutRow * matrixInColM * sizeof (double));
475 matrix_tran (matrixOutDetailPre, matrixInColM, matrixOutRow,
476 matrixOutDetailTemp, matrixInColM, matrixOutRow);
477 free (matrixOutDetailPre);
478 //free (matrixInRowDetailM);
479 //free (matrixInDetailM);
480
481 /* Final Inverse Transform */
482 filterOutLength = 2 * matrixInColM + filterLen - 1;
483 matrixOutPre = malloc (matrixOutRow * matrixOutCol * sizeof (double));
484 for (row = 0; row < matrixOutRow; row++)
485 idwt_neo ((matrixOutApproxTemp + row * matrixInColM),
486 (matrixOutDetailTemp + row * matrixInColM),
487 matrixInColM, lowRe, hiRe, filterLen,
488 (matrixOutPre + row * matrixOutCol),
489 matrixOutCol);
490 free (matrixOutApproxTemp);
491 free (matrixOutDetailTemp);
492 matrix_tran (matrixOutPre, matrixOutRow, matrixOutCol, matrixOut,
493 matrixOutRow, matrixOutCol);
494 free (matrixOutPre);
495
496 return;
497 }
498
499
500 void
wave_mem_cal(int * pLen,int stride,int * total)501 wave_mem_cal (int *pLen, int stride, int *total)
502 {
503 int count = 0;
504 (*total) = 4 * (pLen[2]) * (pLen[3]);
505 for (count = 2; count < (stride+1); count++)
506 (*total) += 3 * (pLen[2 * count]) * (pLen[2 * count + 1]);
507 return;
508 }
509
510 void
matrix_wavedec_len_cal(int matrixInRow,int matrixInCol,int stride,int filterLen,int * pLen)511 matrix_wavedec_len_cal (int matrixInRow, int matrixInCol, int stride,
512 int filterLen, int *pLen)
513 {
514 int row;
515 int filterOutLength;
516
517 pLen[(stride+1) * 2] = matrixInRow;
518 pLen[(stride+1) * 2 + 1] = matrixInCol;
519
520 if (getdwtMode()!=PER)
521 {
522 for (row = stride; row > 0; row--)
523 {
524 filterOutLength = pLen[(row + 1) * 2] + filterLen - 1;
525 pLen[row*2] = filterOutLength/2;
526 filterOutLength = pLen[(row + 1) * 2 + 1] + filterLen - 1;
527 pLen[row*2+1] = filterOutLength/2;
528 }
529 pLen[0] = pLen[2];
530 pLen[1] = pLen[3];
531 }
532 else
533 {
534 for (row = stride; row > 0; row--)
535 {
536 //filterOutLength = pLen[(row + 1) * 2] + filterLen - 1;
537 pLen[row*2] = (int)ceil(((double)(pLen[(row + 1) * 2])) / 2.0);
538 //filterOutLength = pLen[(row + 1) * 2 + 1] + filterLen - 1;
539 pLen[row*2+1] = (int)ceil(((double)(pLen[(row + 1) * 2 + 1])) / 2.0);
540 }
541 pLen[0] = pLen[2];
542 pLen[1] = pLen[3];
543
544 }
545 return;
546 }
547
548 void
matrix_locate(int stride,int * pLen,int * pH,int * pV,int * pD)549 matrix_locate (int stride, int *pLen, int *pH, int *pV, int *pD)
550 {
551 int count;
552 int area, acre;
553
554 /* *pH = (*pLen)*(*(pLen+1)); */
555 pH[0] = pLen[0] * pLen[1];
556 /* *pV = 2 * (*pH); */
557 pV[0] = 2 * pH[0];
558 /* *pD = 3 * (*pH); */
559 pD[0] = 3 * pH[0];
560
561 for (count = 1; count < stride; count++)
562 {
563 /* area = (*(pLen+(count-1)*2)) * (*(pLen+(count-1)*2+1)); */
564 area = pLen[2 * count] * pLen[2 * count + 1];
565 /* acre = (*(pLen+count*2)) * (*(pLen+count*2+1)); */
566 acre = pLen[2 * (count + 1)] * pLen[2 * (count + 1) + 1];
567 /* *(pH+count) = *(pH+count-1) + 2 * area; */
568 pH[count] = pH[count - 1] + 3 * area;
569 /* *(pV+count) = *(pV+count-1) + area + acre; */
570 pV[count] = pV[count - 1] + 2 * area + acre;
571 /* *(pD+count) = *(pD+count-1) + 2 * acre; */
572 pD[count] = pD[count - 1] + area + 2 * acre;
573 }
574 /* *total = *(pD+stride-1)+ acre - 1; */
575
576 return;
577 }
578
579
580 void
wavedec2(double * matrixIn,int matrixInRow,int matrixInCol,double * lowDe,double * hiDe,int filterLen,int * pLen,double * coef,int sigOutLength,int stride,extend_method extMethod)581 wavedec2 (double *matrixIn, int matrixInRow, int matrixInCol,
582 double *lowDe, double *hiDe, int filterLen,
583 int *pLen, double *coef, int sigOutLength, int stride,
584 extend_method extMethod)
585 {
586 int count, row, col;
587 int *pH, *pV, *pD;
588 double *matrixInTemp;
589 double *matrixOutApproxTemp;
590 //printf("waverec2 ok1");
591 matrixInTemp = malloc ((pLen[(stride+1) * 2]) *(pLen[(stride+1) * 2 + 1]) * sizeof (double));
592 matrixOutApproxTemp = malloc ((pLen[stride * 2]) *(pLen[stride * 2 + 1]) * sizeof (double));
593 pH = malloc (stride * sizeof (int));
594 pV = malloc (stride * sizeof (int));
595 pD = malloc (stride * sizeof (int));
596 matrix_locate (stride, pLen, pH, pV, pD);
597
598 for (row = 0; row < pLen[(stride+1) * 2]; row++)
599 {
600 for (col = 0; col < pLen[(stride+1) * 2 + 1]; col++)
601 matrixInTemp[col + row * (pLen[(stride+1) * 2 + 1])] =
602 matrixIn[col + row * (pLen[(stride+1) * 2 + 1])];
603 }
604 //printf("waverec2 ok2");
605 for (count = stride - 1; count >= 0; count--)
606 {
607 dwt2D_neo (matrixInTemp, pLen[2 * (count + 2)],
608 pLen[2 * (count + 2) + 1], matrixOutApproxTemp,
609 (pH[count] + coef), (pV[count] + coef),
610 (pD[count] + coef), pLen[2 * (count + 1)],
611 pLen[2 * (count + 1) + 1], lowDe,
612 hiDe, filterLen, extMethod);
613 //swtdwt2Dex (matrixInTemp, pLen[2 * (count + 1)],
614 // pLen[2 * (count + 1) + 1], matrixOutApproxTemp,
615 // (pH[count] + coef), (pV[count] + coef),
616 // (pD[count] + coef), waveType, mode);
617 for (row = 0; row < pLen[2 * (count+1)]; row++)
618 {
619 for (col = 0; col < pLen[2 * (count+1) + 1]; col++)
620 matrixInTemp[col + row * (pLen[2 * (count+1) + 1])] =
621 matrixOutApproxTemp[col + row * (pLen[2 * (count+1) + 1])];
622 // printf("");
623 }
624 }
625 free (matrixInTemp);
626 free (pH);
627 free (pV);
628 free (pD);
629 // printf("waverec2 ok3");
630 for (row = 0; row < pLen[0]; row++)
631 {
632 for (col = 0; col < pLen[1]; col++)
633 coef[col + row * (pLen[1])] =
634 matrixOutApproxTemp[col + row * (pLen[1])];
635 }
636 free (matrixOutApproxTemp);
637 return;
638 }
639
640 void
waverec2(double * coef,int sigInLength,double * lowRe,double * hiRe,int filterLen,double * matrixOut,int matrixOutRow,int matrixOutCol,int * pLen,int stride,extend_method extMethod)641 waverec2 (double *coef, int sigInLength, double *lowRe, double *hiRe,
642 int filterLen, double *matrixOut, int matrixOutRow,
643 int matrixOutCol, int *pLen, int stride,
644 extend_method extMethod)
645 {
646 int count, row, col;
647 int *pH, *pV, *pD;
648 double *matrixOutTemp;
649 double *matrixInApproxTemp;
650
651 matrixOutTemp = malloc ((pLen[(stride+1) * 2]) *
652 (pLen[(stride+1) * 2 + 1]) * sizeof (double));
653 matrixInApproxTemp = malloc ((pLen[(stride+1) * 2]) *
654 (pLen[(stride+1) * 2 + 1]) * sizeof (double));
655
656 pH = malloc (stride * sizeof (int));
657 pV = malloc (stride * sizeof (int));
658 pD = malloc (stride * sizeof (int));
659 matrix_locate (stride, pLen, pH, pV, pD);
660
661 for (row = 0; row < pLen[0]; row++)
662 {
663 for (col = 0; col < pLen[1]; col++)
664 matrixInApproxTemp[col + row * (pLen[1])] =
665 coef[col + row * (pLen[1])];
666 }
667
668 for (count = 0; count < stride; count++)
669 {
670 //idwt2D (matrixInApproxTemp, (coef + pH[count]),
671 // (coef + pV[count]), (coef + pD[count]),
672 //pLen[(count+1) * 2], pLen[(count+1) * 2 + 1],
673 //lowRe, hiRe, filterLen, matrixOutTemp,
674 //pLen[(count + 2) * 2], pLen[(count + 2) * 2 + 1],
675 //extMethod);
676 idwt2D_neo (matrixInApproxTemp, (coef + pH[count]),
677 (coef + pV[count]), (coef + pD[count]),
678 pLen[(count+1) * 2], pLen[(count+1) * 2 + 1],
679 lowRe, hiRe, filterLen, matrixOutTemp,
680 pLen[(count + 2) * 2], pLen[(count + 2) * 2 + 1]);
681 //swtidwt2Dex (matrixInApproxTemp, (coef + pH[count]),
682 // (coef + pV[count]), (coef + pD[count]),
683 // pLen[count * 2], pLen[count * 2 + 1], matrixOutTemp,
684 // pLen[(count + 1) * 2], pLen[(count + 1) * 2 + 1], waveType,
685 // mode);
686 for (row = 0; row < pLen[(count + 2) * 2]; row++)
687 {
688 for (col = 0; col < pLen[(count + 2) * 2 + 1]; col++)
689 matrixInApproxTemp[col + row * (pLen[(count + 2) * 2 + 1])] =
690 matrixOutTemp[col + row * (pLen[(count + 2) * 2 + 1])];
691 }
692 }
693
694 for (row = 0; row < pLen[(stride+1) * 2]; row++)
695 {
696 for (col = 0; col < pLen[(stride+1) * 2 + 1]; col++)
697 matrixOut[col + row * (pLen[(stride+1) * 2 + 1])] =
698 matrixInApproxTemp[col + row * (pLen[(stride+1) * 2 + 1])];
699 }
700
701
702 free (pH);
703 free (pV);
704 free (pD);
705 free (matrixInApproxTemp);
706 free (matrixOutTemp);
707 return;
708 }
709
710 void
wenergy_2output(double * coef,int sigInLength,int * pLen,double * ae,double * de,int deLength,int stride)711 wenergy_2output (double *coef, int sigInLength, int *pLen,
712 double *ae, double *de, int deLength, int stride)
713 {
714 int count, countt, startp, endp;
715 int *pH, *pV, *pD;
716 double ene;
717
718 ene = 0;
719 for(count=0;count<sigInLength;count++)
720 ene += coef[count]*coef[count];
721
722 *ae = 0;
723 for(count=0;count<pLen[0]*pLen[1];count++)
724 *ae += coef[count]*coef[count];
725 *ae = (*ae)*100/ene;
726
727 pH = malloc (stride * sizeof (int));
728 pV = malloc (stride * sizeof (int));
729 pD = malloc (stride * sizeof (int));
730 matrix_locate (stride, pLen, pH, pV, pD);
731
732 for(count=0;count<stride;count++)
733 {
734 startp = pH[count];
735 endp = pH[count] + 3 * pLen[(count+1)*2] * pLen[(count+1)*2+1];
736 de[count] = 0;
737 for(countt=startp;countt<endp;countt++)
738 de[count] += coef[countt]*coef[countt];
739 de[count] = de[count]*100/ene;
740 }
741 free(pH);
742 free(pV);
743 free(pD);
744 return;
745 }
746
747
748 void
wenergy_4output(double * coef,int sigInLength,int * pLen,double * ae,double * he,double * ve,double * de,int deLength,int stride)749 wenergy_4output (double *coef, int sigInLength, int *pLen,
750 double *ae, double *he, double *ve, double *de,
751 int deLength, int stride)
752 {
753 int count, countt, *pH, *pV, *pD;
754 int starth, startv, startd, endh, endv, endd;
755 double ene;
756 ene = 0;
757 for(count=0;count<sigInLength;count++)
758 ene += coef[count]*coef[count];
759
760 *ae = 0;
761 for(count=0;count<pLen[0]*pLen[1];count++)
762 *ae += coef[count]*coef[count];
763 *ae = (*ae)*100/ene;
764
765 pH = malloc (stride * sizeof (int));
766 pV = malloc (stride * sizeof (int));
767 pD = malloc (stride * sizeof (int));
768 matrix_locate (stride, pLen, pH, pV, pD);
769
770 for(count=0;count<stride;count++)
771 {
772 /* horizontal */
773 starth = pH[count];
774 endh = pH[count] + pLen[(count+1)*2] * pLen[(count+1)*2+1];
775 he[count] = 0;
776 for(countt=starth;countt<endh;countt++)
777 he[count] += coef[countt]*coef[countt];
778 he[count] = he[count]*100/ene;
779 /* vertical */
780 startv = pV[count];
781 endv = pV[count] + pLen[(count+1)*2] * pLen[(count+1)*2+1];
782 ve[count] = 0;
783 for(countt=startv;countt<endv;countt++)
784 ve[count] += coef[countt]*coef[countt];
785 ve[count] = ve[count]*100/ene;
786 /* detail */
787 startd = pD[count];
788 endd = pD[count] + pLen[(count+1)*2] * pLen[(count+1)*2+1];
789 de[count] = 0;
790 for(countt=startd;countt<endd;countt++)
791 de[count] += coef[countt]*coef[countt];
792 de[count] = de[count]*100/ene;
793 }
794 free(pH);
795 free(pV);
796 free(pD);
797
798 return;
799 }
800
801
802 void
detcoef2(double * coef,int sigInLength,double * coefOut,int sigOutLength,int * pLen,int stride,int level,char * coefType)803 detcoef2 (double *coef, int sigInLength, double *coefOut,
804 int sigOutLength, int *pLen, int stride, int level,
805 char *coefType)
806 {
807 int row, col, sta;
808 int *pH, *pV, *pD;
809
810 pH = malloc (stride * sizeof (int));
811 pV = malloc (stride * sizeof (int));
812 pD = malloc (stride * sizeof (int));
813 matrix_locate (stride, pLen, pH, pV, pD);
814 if (strcmp (coefType, "h") == 0)
815 {
816 sta = pH[stride - level];
817 }
818 if (strcmp (coefType, "v") == 0)
819 {
820 sta = pV[stride - level];
821 }
822 if (strcmp (coefType, "d") == 0)
823 {
824 sta = pD[stride - level];
825 }
826
827 for (row = 0; row < pLen[(stride - level + 1) * 2]; row++)
828 {
829 for (col = 0; col < pLen[(stride - level + 1) * 2 + 1]; col++)
830 coefOut[col + row * (pLen[(stride - level + 1) * 2 + 1])] =
831 coef[sta + col + row * (pLen[(stride - level + 1) * 2 + 1])];
832 }
833 free (pH);
834 free (pV);
835 free (pD);
836
837 return;
838 }
839
840
841 void
appcoef2(double * coef,int sigInLength,double * lowRe,double * hiRe,int filterLen,double * coefOut,int matrixOutRow,int matrixOutCol,int * pLen,int stride,int level,extend_method extMethod)842 appcoef2 (double *coef, int sigInLength, double *lowRe,
843 double *hiRe, int filterLen, double *coefOut,
844 int matrixOutRow, int matrixOutCol, int *pLen,
845 int stride, int level, extend_method extMethod)
846 {
847 int count;
848
849 if (level == stride)
850 {
851 for (count = 0; count < ((pLen[0]) * (pLen[1])); count++)
852 coefOut[count] = coef[count];
853 }
854 else
855 {
856 waverec2 (coef, sigInLength, lowRe, hiRe, filterLen,
857 coefOut, matrixOutRow, matrixOutCol,
858 pLen, stride-level, extMethod);
859 }
860
861 return;
862 }
863
864
865 void
wrcoef2(double * coef,int sigInLength,double * lowRe,double * hiRe,int filterLen,double * matrixOut,int matrixOutRow,int matrixOutCol,int * pLen,int stride,int level,char * type,extend_method extMethod)866 wrcoef2 (double *coef, int sigInLength, double *lowRe,
867 double *hiRe, int filterLen, double *matrixOut,
868 int matrixOutRow, int matrixOutCol, int *pLen,
869 int stride, int level, char *type, extend_method extMethod)
870 {
871 int count, total, sta, si;
872 double *coefTemp;
873 int *pH, *pV, *pD;
874
875 wave_mem_cal (pLen, stride, &total);
876 coefTemp = malloc (total * sizeof (double));
877
878 pH = malloc (stride * sizeof (int));
879 pV = malloc (stride * sizeof (int));
880 pD = malloc (stride * sizeof (int));
881 matrix_locate (stride, pLen, pH, pV, pD);
882
883 for (count = 0; count < total; count++)
884 coefTemp[count] = 0;
885
886 if (strcmp (type, "h") == 0)
887 {
888 sta = pH[stride - level];
889 si = (pLen[(stride - level + 1) * 2]) * (pLen[(stride - level + 1) * 2 + 1]);
890 }
891 if (strcmp (type, "v") == 0)
892 {
893 sta = pV[stride - level];
894 si = (pLen[(stride - level + 1) * 2]) * (pLen[(stride - level + 1) * 2 + 1]);
895 }
896 if (strcmp (type, "d") == 0)
897 {
898 sta = pD[stride - level];
899 si = (pLen[(stride - level + 1) * 2]) * (pLen[(stride - level + 1) * 2 + 1]);
900 }
901 if (strcmp (type, "a") == 0)
902 {
903 sta = 0;
904 si = (pLen[0]) * (pLen[1]);
905 if (level != stride)
906 {
907 for (count = 1; count <= (stride - level); count++)
908 si += 3 * (pLen[(count) * 2]) * (pLen[(count) * 2 + 1]);
909 }
910 }
911
912 for (count = sta; count < (si + sta); count++)
913 coefTemp[count] = coef[count];
914
915 waverec2 (coefTemp, sigInLength, lowRe, hiRe, filterLen,
916 matrixOut, matrixOutRow, matrixOutCol,
917 pLen, stride, extMethod);
918
919 //waverec2 (coefTemp, pLen, stride, waveType, matrixOut, mode);
920
921 free (pH);
922 free (pV);
923 free (pD);
924 free (coefTemp);
925 return;
926 }
927
928
929 void
upwlev2(double * coef,int sigInLength,double * lowRe,double * hiRe,int filterLen,int * pLen,int matrixRow,int matrixCol,double * approx,int approxLen,double * newCoef,int newCoefLen,int * newLenMatrix,int lenMatrixRow,int lenMatrixCol,int stride,extend_method extMethod)930 upwlev2 (double *coef, int sigInLength, double *lowRe, double *hiRe,
931 int filterLen, int *pLen, int matrixRow, int matrixCol,
932 double *approx, int approxLen, double *newCoef,
933 int newCoefLen, int *newLenMatrix, int lenMatrixRow,
934 int lenMatrixCol, int stride, extend_method extMethod)
935 {
936 int count, *le, *pH, *pV, *pD, row, col, pos;
937
938 for(count=0;count<approxLen;count++)
939 approx[count]=coef[count];
940
941 le = malloc((matrixRow-1)*matrixCol*sizeof(int));
942 for(count=stride+1;count>1;count--)
943 {
944 le[(count-1)*2] = pLen[count*2];
945 le[(count-1)*2+1] = pLen[count*2+1];
946 }
947 le[0] = pLen[4];
948 le[1] = pLen[5];
949
950 //for(count=0;count<=stride;count++)
951 //sciprint("%d %d\n",le[count*2],le[count*2+1]);
952 for (col = 0; col < matrixCol; col++)
953 {
954 for (row = 0; row < (matrixRow-1); row++)
955 {
956 newLenMatrix[row + col * (matrixRow-1)] =
957 le[col + row * matrixCol];
958 }
959 }
960
961 //for (row = 0; row < matrixCol; row++)
962 //{
963 //for (col = 0; col < (matrixRow-1); col++)
964 // newLenMatrix[row + col * matrixCol] = le[col + row * matrixRow];
965 //}
966 //matrix_tran(le,matrixRow-1,matrixCol,newLenMatrix,
967 // matrixRow-1,matrixCol);
968 free(le);
969
970 pH = malloc (stride * sizeof (int));
971 pV = malloc (stride * sizeof (int));
972 pD = malloc (stride * sizeof (int));
973 matrix_locate (stride, pLen, pH, pV, pD);
974
975 pos = 4*pH[0];
976 for(count=sigInLength-1;count>=pos;count--)
977 newCoef[count-sigInLength+newCoefLen] = coef[count];
978
979
980 //idwt2D (coef, coef+pH[0], coef+pV[0], coef+pD[0],
981 // pLen[0], pLen[1], lowRe, hiRe, filterLen, newCoef,
982 //pLen[4], pLen[5], extMethod);
983 idwt2D_neo (coef, coef+pH[0], coef+pV[0], coef+pD[0],
984 pLen[0], pLen[1], lowRe, hiRe, filterLen, newCoef,
985 pLen[4], pLen[5]);
986
987
988 free(pH);
989 free(pV);
990 free(pD);
991 return;
992 }
993
994
995 void
upcoef2(double * matrixIn,int matrixInRow,int matrixInCol,double * lowRe,double * hiRe,int filterLen,double * matrixOut,int matrixOutRow,int matrixOutCol,int matrixOutDefaultRow,int matrixOutDefaultCol,int step,char * type)996 upcoef2 (double *matrixIn, int matrixInRow, int matrixInCol,
997 double *lowRe, double *hiRe, int filterLen,
998 double *matrixOut, int matrixOutRow, int matrixOutCol,
999 int matrixOutDefaultRow, int matrixOutDefaultCol,
1000 int step, char *type)//, extend_method extMethod)
1001 {
1002 double *vo, *matrixOutTemp, *matrixOutPre;
1003 int matrixOutTempRow, matrixOutTempCol, rowLeng, colLeng,count, count1;
1004
1005 // matrixOutTempRow = 2*matrixInRow - filterLen + 2;
1006 // matrixOutTempCol = 2*matrixInCol - filterLen + 2;
1007 matrixOutTempRow = 2*matrixInRow + filterLen - 2;
1008 matrixOutTempCol = 2*matrixInCol + filterLen - 2;
1009
1010 vo = malloc(matrixInRow*matrixInCol*sizeof(double));
1011 for(count = 0;count<matrixInRow*matrixInCol;count++)
1012 vo[count] = 0;
1013 matrixOutTemp = malloc(matrixOutDefaultRow*
1014 matrixOutDefaultCol*sizeof(double));
1015 for(count=0;count<matrixOutDefaultRow*matrixOutDefaultCol;count++)
1016 matrixOutTemp[count] = 0;
1017
1018 if (!strcmp(type,"a"))
1019 {
1020 //idwt2D (matrixIn, vo, vo, vo, matrixInRow, matrixInCol,
1021 // lowRe, hiRe, filterLen, matrixOutTemp,
1022 //matrixOutTempRow, matrixOutTempCol, extMethod);
1023 idwt2D_neo (matrixIn, vo, vo, vo, matrixInRow, matrixInCol,
1024 lowRe, hiRe, filterLen, matrixOutTemp,
1025 matrixOutTempRow, matrixOutTempCol);
1026 }
1027
1028 if (!strcmp(type,"h"))
1029 {
1030 //idwt2D (vo, matrixIn, vo, vo, matrixInRow, matrixInCol,
1031 // lowRe, hiRe, filterLen, matrixOutTemp,
1032 //matrixOutTempRow, matrixOutTempCol, extMethod);
1033 idwt2D_neo (vo, matrixIn, vo, vo, matrixInRow, matrixInCol,
1034 lowRe, hiRe, filterLen, matrixOutTemp,
1035 matrixOutTempRow, matrixOutTempCol);
1036 }
1037
1038 if (!strcmp(type,"v"))
1039 {
1040 //idwt2D (vo, vo, matrixIn, vo, matrixInRow, matrixInCol,
1041 // lowRe, hiRe, filterLen, matrixOutTemp,
1042 //matrixOutTempRow, matrixOutTempCol, extMethod);
1043 idwt2D_neo (vo, vo, matrixIn, vo, matrixInRow, matrixInCol,
1044 lowRe, hiRe, filterLen, matrixOutTemp,
1045 matrixOutTempRow, matrixOutTempCol);
1046 }
1047
1048 if (!strcmp(type,"d"))
1049 {
1050 //idwt2D (vo, vo, vo, matrixIn, matrixInRow, matrixInCol,
1051 // lowRe, hiRe, filterLen, matrixOutTemp,
1052 //matrixOutTempRow, matrixOutTempCol, extMethod);
1053 idwt2D_neo (vo, vo, vo, matrixIn, matrixInRow, matrixInCol,
1054 lowRe, hiRe, filterLen, matrixOutTemp,
1055 matrixOutTempRow, matrixOutTempCol);
1056 }
1057
1058 free(vo);
1059
1060 if (step>1)
1061 {
1062 matrixOutPre = malloc(matrixOutDefaultRow*
1063 matrixOutDefaultCol*sizeof(double));
1064 for(count=0;count<matrixOutDefaultRow*
1065 matrixOutDefaultCol;count++)
1066 matrixOutPre[count] = 0;
1067 rowLeng = matrixOutTempRow;
1068 colLeng = matrixOutTempCol;
1069 for(count=0;count<step-1;count++)
1070 {
1071 //idwt2D (matrixOutTemp, vo, vo, vo, rowLeng, colLeng,
1072 // lowRe, hiRe, filterLen, matrixOutPre,
1073 //rowLeng*2-filterLen+2, colLeng*2-filterLen+2,
1074 //extMethod);
1075 vo = malloc(rowLeng*colLeng*sizeof(double));
1076 for(count1 = 0;count1<rowLeng*colLeng;count1++)
1077 vo[count1] = 0;
1078 idwt2D_neo (matrixOutTemp, vo, vo, vo, rowLeng, colLeng,
1079 lowRe, hiRe, filterLen, matrixOutPre,
1080 rowLeng*2+filterLen-2, colLeng*2+filterLen-2);
1081 rowLeng = rowLeng*2 + filterLen - 2;
1082 colLeng = colLeng*2 + filterLen - 2;
1083 // idwt2D_neo (matrixOutTemp, vo, vo, vo, rowLeng, colLeng,
1084 // lowRe, hiRe, filterLen, matrixOutPre,
1085 // rowLeng*2-filterLen+2, colLeng*2-filterLen+2);
1086 // rowLeng = rowLeng*2 - filterLen + 2;
1087 // colLeng = colLeng*2 - filterLen + 2;
1088 verbatim_copy (matrixOutPre, rowLeng*colLeng,
1089 matrixOutTemp, rowLeng*colLeng);
1090 free(vo);
1091 }
1092 matrixOutTempRow = rowLeng;
1093 matrixOutTempCol = colLeng;
1094 free(matrixOutPre);
1095 }
1096
1097 wkeep_2D_center (matrixOutTemp, matrixOutDefaultRow,
1098 matrixOutDefaultCol, matrixOut,
1099 matrixOutRow, matrixOutCol);
1100
1101
1102 //free(vo);
1103 free(matrixOutTemp);
1104
1105 return;
1106 }
1107
1108 void
wavedec2a(double * matrixIn,int matrixInRow,int matrixInCol,double * lowDeR,double * hiDeR,double * lowDeC,double * hiDeC,int filterLen,int * pLen,double * coef,int sigOutLength,int stride,extend_method extMethod)1109 wavedec2a (double *matrixIn, int matrixInRow, int matrixInCol,
1110 double *lowDeR, double *hiDeR, double *lowDeC,
1111 double *hiDeC, int filterLen, int *pLen,
1112 double *coef, int sigOutLength, int stride,
1113 extend_method extMethod)
1114 {
1115 int count, row, col;
1116 int *pH, *pV, *pD;
1117 double *matrixInTemp;
1118 double *matrixOutApproxTemp;
1119
1120 matrixInTemp = malloc ((pLen[(stride+1) * 2]) *
1121 (pLen[(stride+1) * 2 + 1]) * sizeof (double));
1122 matrixOutApproxTemp = malloc ((pLen[stride * 2]) *
1123 (pLen[stride * 2 + 1]) *
1124 sizeof (double));
1125 pH = malloc (stride * sizeof (int));
1126 pV = malloc (stride * sizeof (int));
1127 pD = malloc (stride * sizeof (int));
1128 matrix_locate (stride, pLen, pH, pV, pD);
1129
1130 for (row = 0; row < pLen[(stride+1) * 2]; row++)
1131 {
1132 for (col = 0; col < pLen[(stride+1) * 2 + 1]; col++)
1133 matrixInTemp[col + row * (pLen[(stride+1) * 2 + 1])] =
1134 matrixIn[col + row * (pLen[(stride+1) * 2 + 1])];
1135 }
1136
1137 for (count = stride - 1; count >= 0; count--)
1138 {
1139 dwt2D_neo_a (matrixInTemp, pLen[2 * (count + 2)],
1140 pLen[2 * (count + 2) + 1], matrixOutApproxTemp,
1141 (pH[count] + coef), (pV[count] + coef),
1142 (pD[count] + coef), pLen[2 * (count + 1)],
1143 pLen[2 * (count + 1) + 1], lowDeR,
1144 hiDeR, lowDeC, hiDeC, filterLen, extMethod);
1145 for (row = 0; row < pLen[2 * (count+1)]; row++)
1146 {
1147 for (col = 0; col < pLen[2 * (count+1) + 1]; col++)
1148 matrixInTemp[col + row * (pLen[2 * (count+1) + 1])] =
1149 matrixOutApproxTemp[col + row * (pLen[2 * (count+1) + 1])];
1150 }
1151 }
1152 free (matrixInTemp);
1153 free (pH);
1154 free (pV);
1155 free (pD);
1156 for (row = 0; row < pLen[0]; row++)
1157 {
1158 for (col = 0; col < pLen[1]; col++)
1159 coef[col + row * (pLen[1])] =
1160 matrixOutApproxTemp[col + row * (pLen[1])];
1161 }
1162 free (matrixOutApproxTemp);
1163 return;
1164 }
1165
1166 void
dwt2D_neo_a(double * matrixIn,int matrixInRow,int matrixInCol,double * matrixOutApprox,double * matrixOutColDetail,double * matrixOutRowDetail,double * matrixOutDetail,int matrixOutRow,int matrixOutCol,double * lowDeR,double * hiDeR,double * lowDeC,double * hiDeC,int filterLen,extend_method extMethod)1167 dwt2D_neo_a (double *matrixIn, int matrixInRow, int matrixInCol,
1168 double *matrixOutApprox, double *matrixOutColDetail,
1169 double *matrixOutRowDetail, double *matrixOutDetail,
1170 int matrixOutRow, int matrixOutCol, double *lowDeR,
1171 double *hiDeR, double *lowDeC, double *hiDeC,
1172 int filterLen, extend_method extMethod)
1173 {
1174 int row, col;
1175 int filterOutLength;
1176 int matrixOutRowM, matrixOutColM, matrixInRowM, matrixInColM;
1177 int matrixOutRowT, matrixOutColT;
1178 double *matrixInPre, *matrixInRo;
1179 double *pMatrixOutApproxPre, *pMatrixOutDetailPre;
1180 double *matrixOutApproxPre, *matrixOutDetailPre;
1181 double *matrixOutApproxM, *matrixOutColDetailM;
1182 double *matrixOutRowDetailM, *matrixOutDetailM;
1183 double *matrixOutApproxT, *matrixOutColDetailT;
1184 double *matrixOutRowDetailT, *matrixOutDetailT;
1185 char c='b';
1186
1187 /* extension */
1188 matrixInRowM = matrixInRow + 2 * filterLen;
1189 matrixInColM = matrixInCol + 2 * filterLen;
1190 if ((extMethod==PER)&&(matrixInRow%2!=0))
1191 matrixInRowM = matrixInRow + 2 * filterLen + 1;
1192 if ((extMethod==PER)&&(matrixInCol%2!=0))
1193 matrixInColM = matrixInCol + 2 * filterLen + 1;
1194
1195 matrixInRo = malloc (matrixInRowM * matrixInColM * sizeof (double));
1196 matrixInPre = malloc(matrixInRowM * matrixInColM * sizeof (double));
1197 wextend_2D (matrixIn, matrixInRow, matrixInCol, matrixInRo,
1198 matrixInRowM, matrixInColM, extMethod, &c, &c);
1199 matrix_tran (matrixInRo, matrixInColM, matrixInRowM,
1200 matrixInPre, matrixInRowM, matrixInColM);
1201 free (matrixInRo);
1202
1203 /* Row DWT */
1204
1205 filterOutLength = matrixInColM + filterLen - 1;
1206 matrixOutColM = filterOutLength;
1207 filterOutLength = matrixInRowM + filterLen - 1;
1208 matrixOutRowM = filterOutLength;
1209
1210 matrixOutApproxPre = malloc (matrixInRowM *
1211 matrixOutColM * sizeof (double));
1212 matrixOutDetailPre = malloc (matrixInRowM *
1213 matrixOutColM * sizeof (double));
1214 filterOutLength = matrixInColM + filterLen - 1;
1215 for (row = 0; row < matrixInRowM; row++)
1216 {
1217 dwt_conv ((matrixInPre + row * matrixInColM), matrixInColM, lowDeR,
1218 hiDeR, filterLen,
1219 (matrixOutApproxPre + row * matrixOutColM),
1220 (matrixOutDetailPre + row * matrixOutColM),
1221 matrixOutColM);
1222 }
1223
1224 free (matrixInPre);
1225
1226 /* transpose */
1227 pMatrixOutApproxPre = malloc (matrixInRowM *
1228 matrixOutColM * sizeof (double));
1229 matrix_tran (matrixOutApproxPre, matrixInRowM, matrixOutColM,
1230 pMatrixOutApproxPre, matrixInRowM, matrixOutColM);
1231 free (matrixOutApproxPre);
1232
1233 pMatrixOutDetailPre = malloc (matrixInRowM *
1234 matrixOutColM * sizeof (double));
1235 matrix_tran (matrixOutDetailPre, matrixInRowM, matrixOutColM,
1236 pMatrixOutDetailPre, matrixInRowM, matrixOutColM);
1237
1238 free (matrixOutDetailPre);
1239
1240 /* Col DWT */
1241 filterOutLength = matrixInRowM + filterLen - 1;
1242 matrixOutApproxM = malloc (matrixOutRowM * matrixOutColM
1243 * sizeof (double));
1244 matrixOutColDetailM = malloc (matrixOutRowM * matrixOutColM
1245 * sizeof (double));
1246
1247 for (col = 0; col < matrixOutColM; col++)
1248 {
1249 dwt_conv ((pMatrixOutApproxPre + col * matrixInRowM), matrixInRowM,
1250 lowDeC, hiDeC, filterLen,
1251 (matrixOutApproxM + col * matrixOutRowM),
1252 (matrixOutColDetailM + col * matrixOutRowM),
1253 matrixOutRowM);//, extMethod);
1254 }
1255 free (pMatrixOutApproxPre);
1256
1257 matrixOutRowT = matrixInRow + filterLen - 1;
1258 matrixOutColT = matrixInCol + filterLen - 1;
1259 if ((extMethod==PER)&&(matrixInRow%2!=0))
1260 matrixOutRowT = matrixInRow + 1;
1261 if ((extMethod==PER)&&(matrixInCol%2!=0))
1262 matrixOutColT = matrixInCol + 1;
1263 if ((extMethod==PER)&&(matrixInRow%2==0))
1264 matrixOutRowT = matrixInRow;
1265 if ((extMethod==PER)&&(matrixInCol%2==0))
1266 matrixOutColT = matrixInCol;
1267 matrixOutApproxT = malloc (matrixOutRowT*matrixOutColT*sizeof(double));
1268 matrixOutColDetailT = malloc (matrixOutRowT*matrixOutColT*sizeof(double));
1269
1270 wkeep_2D_center (matrixOutApproxM, matrixOutRowM, matrixOutColM,
1271 matrixOutApproxT, matrixOutRowT, matrixOutColT);
1272 free (matrixOutApproxM);
1273 dyaddown_2D_keep_even (matrixOutApproxT, matrixOutRowT,
1274 matrixOutColT, matrixOutApprox,
1275 matrixOutRow, matrixOutCol);
1276 free(matrixOutApproxT);
1277
1278 wkeep_2D_center (matrixOutColDetailM, matrixOutRowM, matrixOutColM,
1279 matrixOutColDetailT, matrixOutRowT, matrixOutColT);
1280 free (matrixOutColDetailM);
1281 dyaddown_2D_keep_even (matrixOutColDetailT, matrixOutRowT,
1282 matrixOutColT, matrixOutColDetail,
1283 matrixOutRow, matrixOutCol);
1284 free(matrixOutColDetailT);
1285
1286
1287
1288
1289 filterOutLength = matrixInRowM + filterLen - 1;
1290 matrixOutRowDetailM = malloc (matrixOutRowM * matrixOutColM
1291 * sizeof (double));
1292 matrixOutDetailM = malloc (matrixOutRowM * matrixOutColM * sizeof (double));
1293 for (col = 0; col < matrixOutColM; col++)
1294 {
1295 dwt_conv ((pMatrixOutDetailPre + col * matrixInRowM), matrixInRowM,
1296 lowDeC, hiDeC, filterLen,
1297 (matrixOutRowDetailM + col * matrixOutRowM),
1298 (matrixOutDetailM + col * matrixOutRowM),
1299 matrixOutRowM);
1300 }
1301 free (pMatrixOutDetailPre);
1302 matrixOutRowDetailT = malloc (matrixOutRowT*matrixOutColT*sizeof(double));
1303 matrixOutDetailT = malloc (matrixOutRowT*matrixOutColT*sizeof(double));
1304
1305
1306 wkeep_2D_center (matrixOutRowDetailM, matrixOutRowM, matrixOutColM,
1307 matrixOutRowDetailT, matrixOutRowT, matrixOutColT);
1308 free (matrixOutRowDetailM);
1309 dyaddown_2D_keep_even (matrixOutRowDetailT, matrixOutRowT,
1310 matrixOutColT, matrixOutRowDetail,
1311 matrixOutRow, matrixOutCol);
1312 free(matrixOutRowDetailT);
1313
1314 wkeep_2D_center (matrixOutDetailM, matrixOutRowM, matrixOutColM,
1315 matrixOutDetailT, matrixOutRowT, matrixOutColT);
1316 free (matrixOutDetailM);
1317 dyaddown_2D_keep_even (matrixOutDetailT, matrixOutRowT,
1318 matrixOutColT, matrixOutDetail,
1319 matrixOutRow, matrixOutCol);
1320 free(matrixOutDetailT);
1321
1322 return;
1323 }
1324
1325 void
waverec2a(double * coef,int sigInLength,double * lowReR,double * hiReR,double * lowReC,double * hiReC,int filterLen,double * matrixOut,int matrixOutRow,int matrixOutCol,int * pLen,int stride,extend_method extMethod)1326 waverec2a (double *coef, int sigInLength, double *lowReR,
1327 double *hiReR, double *lowReC, double *hiReC,
1328 int filterLen, double *matrixOut, int matrixOutRow,
1329 int matrixOutCol, int *pLen, int stride,
1330 extend_method extMethod)
1331 {
1332 int count, row, col;
1333 int *pH, *pV, *pD;
1334 double *matrixOutTemp;
1335 double *matrixInApproxTemp;
1336
1337 matrixOutTemp = malloc ((pLen[(stride+1) * 2]) *
1338 (pLen[(stride+1) * 2 + 1]) * sizeof (double));
1339 matrixInApproxTemp = malloc ((pLen[(stride+1) * 2]) *
1340 (pLen[(stride+1) * 2 + 1]) * sizeof (double));
1341
1342 pH = malloc (stride * sizeof (int));
1343 pV = malloc (stride * sizeof (int));
1344 pD = malloc (stride * sizeof (int));
1345 matrix_locate (stride, pLen, pH, pV, pD);
1346
1347 for (row = 0; row < pLen[0]; row++)
1348 {
1349 for (col = 0; col < pLen[1]; col++)
1350 matrixInApproxTemp[col + row * (pLen[1])] =
1351 coef[col + row * (pLen[1])];
1352 }
1353
1354 for (count = 0; count < stride; count++)
1355 {
1356 idwt2D_neo_a (matrixInApproxTemp, (coef + pH[count]),
1357 (coef + pV[count]), (coef + pD[count]),
1358 pLen[(count+1) * 2], pLen[(count+1) * 2 + 1],
1359 lowReR, hiReR, lowReC, hiReC, filterLen,
1360 matrixOutTemp, pLen[(count + 2) * 2],
1361 pLen[(count + 2) * 2 + 1]);
1362
1363 for (row = 0; row < pLen[(count + 2) * 2]; row++)
1364 {
1365 for (col = 0; col < pLen[(count + 2) * 2 + 1]; col++)
1366 matrixInApproxTemp[col + row * (pLen[(count + 2) * 2 + 1])] =
1367 matrixOutTemp[col + row * (pLen[(count + 2) * 2 + 1])];
1368 }
1369 }
1370
1371 for (row = 0; row < pLen[(stride+1) * 2]; row++)
1372 {
1373 for (col = 0; col < pLen[(stride+1) * 2 + 1]; col++)
1374 matrixOut[col + row * (pLen[(stride+1) * 2 + 1])] =
1375 matrixInApproxTemp[col + row * (pLen[(stride+1) * 2 + 1])];
1376 }
1377
1378
1379 free (pH);
1380 free (pV);
1381 free (pD);
1382 free (matrixInApproxTemp);
1383 free (matrixOutTemp);
1384 return;
1385 }
1386
1387 void
idwt2D_neo_a(double * matrixInApprox,double * matrixInColDetail,double * matrixInRowDetail,double * matrixInDetail,int matrixInRow,int matrixInCol,double * lowReR,double * hiReR,double * lowReC,double * hiReC,int filterLen,double * matrixOut,int matrixOutRow,int matrixOutCol)1388 idwt2D_neo_a (double *matrixInApprox, double *matrixInColDetail,
1389 double *matrixInRowDetail, double *matrixInDetail,
1390 int matrixInRow, int matrixInCol, double *lowReR,
1391 double *hiReR, double *lowReC, double *hiReC,
1392 int filterLen, double *matrixOut,
1393 int matrixOutRow, int matrixOutCol)
1394 {
1395 int row, col;
1396 int filterOutLength;
1397 int matrixInRowM, matrixInColM;
1398 char c='b';
1399 double *matrixOutApproxPre, *matrixOutDetailPre;
1400 double *matrixOutApproxTemp, *matrixOutDetailTemp;
1401 double *matrixOutPre;
1402
1403 matrixInRowM = matrixInRow;// + 2 * (filterLen - 1);
1404 matrixInColM = matrixInCol;// + 2 * (filterLen - 1);
1405
1406
1407 filterOutLength = 2 * matrixInRowM + filterLen - 1;
1408
1409 /* Approx Calculation */
1410 matrixOutApproxPre = malloc (matrixOutRow * matrixInColM * sizeof (double));
1411 matrixOutApproxTemp =
1412 malloc (matrixOutRow * matrixInColM * sizeof (double));
1413 for (col = 0; col < matrixInColM; col++)
1414 idwt_neo ((matrixInApprox + col * matrixInRowM),
1415 (matrixInColDetail + col * matrixInRowM),
1416 matrixInRowM, lowReC, hiReC, filterLen,
1417 (matrixOutApproxPre + col * matrixOutRow),
1418 matrixOutRow);
1419 matrix_tran (matrixOutApproxPre, matrixInColM, matrixOutRow,
1420 matrixOutApproxTemp, matrixInColM, matrixOutRow);
1421 free (matrixOutApproxPre);
1422
1423 /* Detail Calculation */
1424 matrixOutDetailPre = malloc (matrixOutRow * matrixInColM * sizeof (double));
1425 for (col = 0; col < matrixInColM; col++)
1426 idwt_neo ((matrixInRowDetail + col * matrixInRowM),
1427 (matrixInDetail + col * matrixInRowM),
1428 matrixInRowM, lowReC, hiReC, filterLen,
1429 (matrixOutDetailPre + col * matrixOutRow),
1430 matrixOutRow);
1431 matrixOutDetailTemp =
1432 malloc (matrixOutRow * matrixInColM * sizeof (double));
1433 matrix_tran (matrixOutDetailPre, matrixInColM, matrixOutRow,
1434 matrixOutDetailTemp, matrixInColM, matrixOutRow);
1435 free (matrixOutDetailPre);
1436
1437 /* Final Inverse Transform */
1438 filterOutLength = 2 * matrixInColM + filterLen - 1;
1439 matrixOutPre = malloc (matrixOutRow * matrixOutCol * sizeof (double));
1440 for (row = 0; row < matrixOutRow; row++)
1441 idwt_neo ((matrixOutApproxTemp + row * matrixInColM),
1442 (matrixOutDetailTemp + row * matrixInColM),
1443 matrixInColM, lowReR, hiReR, filterLen,
1444 (matrixOutPre + row * matrixOutCol),
1445 matrixOutCol);
1446 free (matrixOutApproxTemp);
1447 free (matrixOutDetailTemp);
1448 matrix_tran (matrixOutPre, matrixOutRow, matrixOutCol, matrixOut,
1449 matrixOutRow, matrixOutCol);
1450 free (matrixOutPre);
1451
1452 return;
1453 }
1454