From d6a523e091c46ba502329579b01df6ff2a7b3fde Mon Sep 17 00:00:00 2001 From: Antonin Descampe Date: Wed, 21 Sep 2005 13:00:51 +0000 Subject: [PATCH] major change in the dwt-module, thanks to Ive (aka Reiner Wahler): thanks a lot ! See note in dwt.c for more details. --- libopenjpeg/dwt.c | 567 ++++++++++++++++++++++++++-------------------- libopenjpeg/dwt.h | 28 +-- libopenjpeg/tcd.c | 69 +++--- 3 files changed, 363 insertions(+), 301 deletions(-) diff --git a/libopenjpeg/dwt.c b/libopenjpeg/dwt.c index 92d006ad..ef567b00 100644 --- a/libopenjpeg/dwt.c +++ b/libopenjpeg/dwt.c @@ -2,6 +2,7 @@ * Copyright (c) 2001-2002, David Janssens * Copyright (c) 2002-2004, Yannick Verschueren * Copyright (c) 2002-2004, Communications and remote sensing Laboratory, Universite catholique de Louvain, Belgium + * Copyright (c) 2005, Reiner Wahler * All rights reserved. * * Redistribution and use in source and binary forms, with or without @@ -26,16 +27,39 @@ * POSSIBILITY OF SUCH DAMAGE. */ +/* + * NOTE: + * This is a modified version of the openjpeg dwt.c file. + * Average speed improvement compared to the original file (measured on + * my own machine, a P4 running at 3.0 GHz): + * 5x3 wavelets about 2 times faster + * 9x7 wavelets about 3 times faster + * for both, encoding and decoding. + * + * The better performance is caused by doing the 1-dimensional DWT + * within a temporary buffer where the data can be accessed sequential + * for both directions, horizontal and vertical. The 2d vertical DWT was + * the major bottleneck in the former version. + * + * I have also removed the "Add Patrick" part because it is not longer + * needed. + * + * 6/6/2005 + * -Ive (aka Reiner Wahler) + * mail: ive@lilysoft.com + */ + + #include "dwt.h" #include "int.h" #include "fix.h" #include "tcd.h" #include -#include +//#include //#include -#define S(i) a[x*(i)*2] -#define D(i) a[x*(1+(i)*2)] +#define S(i) a[(i)*2] +#define D(i) a[(1+(i)*2)] #define S_(i) ((i)<0?S(0):((i)>=sn?S(sn-1):S(i))) #define D_(i) ((i)<0?D(0):((i)>=dn?D(dn-1):D(i))) /* new */ @@ -62,211 +86,230 @@ double dwt_norms_real[4][10] = { {2.080, 3.865, 8.307, 17.18, 34.71, 69.59, 139.3, 278.6, 557.2} }; -/* Add Patrick */ -static int *b = NULL; -static int lastSizeOfB = 0; -/* */ -/* Cleaning memory. */ -/* */ - -void dwt_clean() -{ - if (b != NULL) { - free(b); - } - b = NULL; - lastSizeOfB = 0; +/* */ +/* Forward lazy transform (horizontal). */ +/* */ +void dwt_deinterleave_h(int *a, int *b, int dn, int sn, int cas) { + int i; + for (i=0; i */ -/* Forward lazy transform. */ -/* */ -void dwt_deinterleave(int *a, int n, int x, int res, int cas) -{ - int dn, sn, i; - sn = res; - dn = n - res; - if (lastSizeOfB != n) { - if (b != NULL) - free(b); - b = (int *) malloc(n * sizeof(int)); - lastSizeOfB = n; - } - - if (cas) { - for (i = 0; i < sn; i++) - b[i] = a[(2 * i + 1) * x]; - for (i = 0; i < dn; i++) - b[sn + i] = a[2 * i * x]; - } else { - for (i = 0; i < sn; i++) - b[i] = a[2 * i * x]; - for (i = 0; i < dn; i++) - b[sn + i] = a[(2 * i + 1) * x]; - } - for (i = 0; i < n; i++) - a[i * x] = b[i]; +/* */ +/* Forward lazy transform (vertical). */ +/* */ +void dwt_deinterleave_v(int *a, int *b, int dn, int sn, int x, int cas) { + int i; + for (i=0; i */ -/* Inverse lazy transform. */ -/* */ -void dwt_interleave(int *a, int n, int x, int res, int cas) -{ - int dn, sn, i; - sn = res; - dn = n - res; - - if (lastSizeOfB != n) { - if (b != NULL) - free(b); - b = (int *) malloc(n * sizeof(int)); - lastSizeOfB = n; - } - - if (cas) { - for (i = 0; i < sn; i++) - b[2 * i + 1] = a[i * x]; - for (i = 0; i < dn; i++) - b[2 * i] = a[(sn + i) * x]; - } else { - for (i = 0; i < sn; i++) - b[2 * i] = a[i * x]; - for (i = 0; i < dn; i++) - b[2 * i + 1] = a[(sn + i) * x]; - } - for (i = 0; i < n; i++) - a[i * x] = b[i]; +/* */ +/* Inverse lazy transform (horizontal). */ +/* */ +void dwt_interleave_h(int *a, int *b, int dn, int sn, int cas) { + int i; +// for (i=0; i */ +/* Inverse lazy transform (vertical). */ +/* */ +void dwt_interleave_v(int *a, int *b, int dn, int sn, int x, int cas) { + int i; +// for (i=0; i */ /* Forward 5-3 wavelet tranform in 1-D. */ /* */ -void dwt_encode_1(int *a, int n, int x, int res, int cas) +void dwt_encode_1(int *a, int dn, int sn, int cas) { - int dn, sn, i = 0; - sn = res; - dn = n - res; + int i; - if (cas) { - if (!sn && dn == 1) /* NEW : CASE ONE ELEMENT */ - S(i) *= 2; - else { - for (i = 0; i < dn; i++) - S(i) -= (DD_(i) + DD_(i - 1)) >> 1; - for (i = 0; i < sn; i++) - D(i) += (SS_(i) + SS_(i + 1) + 2) >> 2; + if (!cas) { + if ((dn > 0) || (sn > 1)) { /* NEW : CASE ONE ELEMENT */ + for (i = 0; i < dn; i++) D(i) -= (S_(i) + S_(i + 1)) >> 1; + for (i = 0; i < sn; i++) S(i) += (D_(i - 1) + D_(i) + 2) >> 2; } } else { - if ((dn > 0) || (sn > 1)) { /* NEW : CASE ONE ELEMENT */ - for (i = 0; i < dn; i++) - D(i) -= (S_(i) + S_(i + 1)) >> 1; - for (i = 0; i < sn; i++) - S(i) += (D_(i - 1) + D_(i) + 2) >> 2; + if (!sn && dn == 1) /* NEW : CASE ONE ELEMENT */ + S(0) *= 2; + else { + for (i = 0; i < dn; i++) S(i) -= (DD_(i) + DD_(i - 1)) >> 1; + for (i = 0; i < sn; i++) D(i) += (SS_(i) + SS_(i + 1) + 2) >> 2; } + } - dwt_deinterleave(a, n, x, res, cas); } /* */ /* Inverse 5-3 wavelet tranform in 1-D. */ -/* */ -void dwt_decode_1(int *a, int n, int x, int res, int cas) +/* */ +void dwt_decode_1(int *a, int dn, int sn, int cas) { - int dn, sn, i = 0; - sn = res; - dn = n - res; + int i; - dwt_interleave(a, n, x, res, cas); - if (cas) { - if (!sn && dn == 1) /* NEW : CASE ONE ELEMENT */ - S(i) /= 2; - else { - for (i = 0; i < sn; i++) - D(i) -= (SS_(i) + SS_(i + 1) + 2) >> 2; - for (i = 0; i < dn; i++) - S(i) += (DD_(i) + DD_(i - 1)) >> 1; + if (!cas) { + if ((dn > 0) || (sn > 1)) { /* NEW : CASE ONE ELEMENT */ + for (i = 0; i < sn; i++) S(i) -= (D_(i - 1) + D_(i) + 2) >> 2; + for (i = 0; i < dn; i++) D(i) += (S_(i) + S_(i + 1)) >> 1; } } else { - if ((dn > 0) || (sn > 1)) { /* NEW : CASE ONE ELEMENT */ - for (i = 0; i < sn; i++) - S(i) -= (D_(i - 1) + D_(i) + 2) >> 2; - for (i = 0; i < dn; i++) - D(i) += (S_(i) + S_(i + 1)) >> 1; + if (!sn && dn == 1) /* NEW : CASE ONE ELEMENT */ + S(0) /= 2; + else { + for (i = 0; i < sn; i++) D(i) -= (SS_(i) + SS_(i + 1) + 2) >> 2; + for (i = 0; i < dn; i++) S(i) += (DD_(i) + DD_(i - 1)) >> 1; } } } + /* */ /* Forward 5-3 wavelet tranform in 2-D. */ /* */ -void dwt_encode(int *a, int w, int h, tcd_tilecomp_t * tilec, int l) +void dwt_encode(tcd_tilecomp_t * tilec) { - int i, j; - int rw; /* width of the resolution level computed */ - int rh; /* heigth of the resolution level computed */ - int rw1; /* width of the resolution level once lower than computed one */ - int rh1; /* height of the resolution level once lower than computed one */ + int i, j, k; + int* a; + int* aj; + int* bj; + int w, l; + + w = tilec->x1-tilec->x0; + l = tilec->numresolutions-1; + a = tilec->data; for (i = 0; i < l; i++) { - int cas_col = 0; /* 0 = non inversion on horizontal filtering 1 = inversion between low-pass and high-pass filtering */ - int cas_row = 0; /* 0 = non inversion on vertical filtering 1 = inversion between low-pass and high-pass filtering */ + int rw; /* width of the resolution level computed */ + int rh; /* heigth of the resolution level computed */ + int rw1; /* width of the resolution level once lower than computed one */ + int rh1; /* height of the resolution level once lower than computed one */ + int cas_col; /* 0 = non inversion on horizontal filtering 1 = inversion between low-pass and high-pass filtering */ + int cas_row; /* 0 = non inversion on vertical filtering 1 = inversion between low-pass and high-pass filtering */ + int dn, sn; + rw = tilec->resolutions[l - i].x1 - tilec->resolutions[l - i].x0; rh = tilec->resolutions[l - i].y1 - tilec->resolutions[l - i].y0; - rw1 = - tilec->resolutions[l - i - 1].x1 - tilec->resolutions[l - i - 1].x0; - rh1 = - tilec->resolutions[l - i - 1].y1 - tilec->resolutions[l - i - 1].y0; + rw1= tilec->resolutions[l - i - 1].x1 - tilec->resolutions[l - i - 1].x0; + rh1= tilec->resolutions[l - i - 1].y1 - tilec->resolutions[l - i - 1].y0; cas_row = tilec->resolutions[l - i].x0 % 2; cas_col = tilec->resolutions[l - i].y0 % 2; - for (j = 0; j < rw; j++) - dwt_encode_1(a + j, rh, w, rh1, cas_col); - for (j = 0; j < rh; j++) - dwt_encode_1(a + j * w, rw, 1, rw1, cas_row); - } - dwt_clean(); + sn = rh1; + dn = rh - rh1; + bj=(int*)malloc(rh*sizeof(int)); + for (j=0; j */ /* Inverse 5-3 wavelet tranform in 2-D. */ /* */ -void dwt_decode(int *a, int w, int h, tcd_tilecomp_t * tilec, int l, - int stop) +void dwt_decode(tcd_tilecomp_t * tilec, int stop) { - int i, j; - int rw; /* width of the resolution level computed */ - int rh; /* heigth of the resolution level computed */ - int rw1; /* width of the resolution level once lower than computed one */ - int rh1; /* height of the resolution level once lower than computed one */ + int i, j, k; + int* a; + int* aj; + int* bj; + int w, l; + + w = tilec->x1-tilec->x0; + l = tilec->numresolutions-1; + a = tilec->data; for (i = l - 1; i >= stop; i--) { - int cas_col = 0; /* 0 = non inversion on horizontal filtering 1 = inversion between low-pass and high-pass filtering */ - int cas_row = 0; /* 0 = non inversion on vertical filtering 1 = inversion between low-pass and high-pass filtering */ + int rw; /* width of the resolution level computed */ + int rh; /* heigth of the resolution level computed */ + int rw1; /* width of the resolution level once lower than computed one */ + int rh1; /* height of the resolution level once lower than computed one */ + int cas_col; /* 0 = non inversion on horizontal filtering 1 = inversion between low-pass and high-pass filtering */ + int cas_row; /* 0 = non inversion on vertical filtering 1 = inversion between low-pass and high-pass filtering */ + int dn, sn; rw = tilec->resolutions[l - i].x1 - tilec->resolutions[l - i].x0; rh = tilec->resolutions[l - i].y1 - tilec->resolutions[l - i].y0; - rw1 = - tilec->resolutions[l - i - 1].x1 - tilec->resolutions[l - i - 1].x0; - rh1 = - tilec->resolutions[l - i - 1].y1 - tilec->resolutions[l - i - 1].y0; + rw1= tilec->resolutions[l - i - 1].x1 - tilec->resolutions[l - i - 1].x0; + rh1= tilec->resolutions[l - i - 1].y1 - tilec->resolutions[l - i - 1].y0; cas_row = tilec->resolutions[l - i].x0 % 2; cas_col = tilec->resolutions[l - i].y0 % 2; - for (j = 0; j < rh; j++) - dwt_decode_1(a + j * w, rw, 1, rw1, cas_row); - for (j = 0; j < rw; j++) - dwt_decode_1(a + j, rh, w, rh1, cas_col); + sn = rw1; + dn = rw - rw1; + bj=(int*)malloc(rw*sizeof(int)); + for (j = 0; j < rh; j++) { + aj = a+j*w; + dwt_interleave_h(aj, bj, dn, sn, cas_row); + dwt_decode_1(bj, dn, sn, cas_row); + for (k = 0; k < rw; k++) aj[k] = bj[k]; + } + free(bj); + + sn = rh1; + dn = rh - rh1; + bj=(int*)malloc(rh*sizeof(int)); + for (j = 0; j < rw; j++) { + aj = a+j; + dwt_interleave_v(aj, bj, dn, sn, w, cas_col); + dwt_decode_1(bj, dn, sn, cas_col); + for (k = 0; k < rh; k++) aj[k * w] = bj[k]; + } + free(bj); + } - dwt_clean(); } + /* */ /* Get gain of 5-3 wavelet transform. */ /* */ @@ -290,84 +333,77 @@ double dwt_getnorm(int level, int orient) /* */ /* Forward 9-7 wavelet transform in 1-D. */ /* */ -void dwt_encode_1_real(int *a, int n, int x, int res, int cas) +void dwt_encode_1_real(int *a, int dn, int sn, int cas) { - int dn, sn, i = 0; - dn = n - res; - sn = res; - - if (cas) { - if ((sn > 0) || (dn > 1)) { /* NEW : CASE ONE ELEMENT */ - for (i = 0; i < dn; i++) - S(i) -= fix_mul(DD_(i) + DD_(i - 1), 12993); - for (i = 0; i < sn; i++) - D(i) -= fix_mul(SS_(i) + SS_(i + 1), 434); - for (i = 0; i < dn; i++) - S(i) += fix_mul(DD_(i) + DD_(i - 1), 7233); - for (i = 0; i < sn; i++) - D(i) += fix_mul(SS_(i) + SS_(i + 1), 3633); - for (i = 0; i < dn; i++) - S(i) = fix_mul(S(i), 5038); /*5038 */ - for (i = 0; i < sn; i++) - D(i) = fix_mul(D(i), 6659); /*6660 */ - } - } else { + int i; + if (!cas) { if ((dn > 0) || (sn > 1)) { /* NEW : CASE ONE ELEMENT */ for (i = 0; i < dn; i++) - D(i) -= fix_mul(S_(i) + S_(i + 1), 12993); + D(i) -= fix_mul(S_(i) + S_(i + 1), 12993); for (i = 0; i < sn; i++) - S(i) -= fix_mul(D_(i - 1) + D_(i), 434); + S(i) -= fix_mul(D_(i - 1) + D_(i), 434); for (i = 0; i < dn; i++) - D(i) += fix_mul(S_(i) + S_(i + 1), 7233); + D(i) += fix_mul(S_(i) + S_(i + 1), 7233); for (i = 0; i < sn; i++) - S(i) += fix_mul(D_(i - 1) + D_(i), 3633); + S(i) += fix_mul(D_(i - 1) + D_(i), 3633); for (i = 0; i < dn; i++) - D(i) = fix_mul(D(i), 5038); /*5038 */ + D(i) = fix_mul(D(i), 5038); /*5038 */ for (i = 0; i < sn; i++) - S(i) = fix_mul(S(i), 6659); /*6660 */ + S(i) = fix_mul(S(i), 6659); /*6660 */ + } + } else { + if ((sn > 0) || (dn > 1)) { /* NEW : CASE ONE ELEMENT */ + for (i = 0; i < dn; i++) + S(i) -= fix_mul(DD_(i) + DD_(i - 1), 12993); + for (i = 0; i < sn; i++) + D(i) -= fix_mul(SS_(i) + SS_(i + 1), 434); + for (i = 0; i < dn; i++) + S(i) += fix_mul(DD_(i) + DD_(i - 1), 7233); + for (i = 0; i < sn; i++) + D(i) += fix_mul(SS_(i) + SS_(i + 1), 3633); + for (i = 0; i < dn; i++) + S(i) = fix_mul(S(i), 5038); /*5038 */ + for (i = 0; i < sn; i++) + D(i) = fix_mul(D(i), 6659); /*6660 */ } } - dwt_deinterleave(a, n, x, res, cas); } /* */ /* Inverse 9-7 wavelet transform in 1-D. */ /* */ -void dwt_decode_1_real(int *a, int n, int x, int res, int cas) +void dwt_decode_1_real(int *a, int dn, int sn, int cas) { - int dn, sn, i = 0; - dn = n - res; - sn = res; - dwt_interleave(a, n, x, res, cas); - if (cas) { - if ((sn > 0) || (dn > 1)) { /* NEW : CASE ONE ELEMENT */ - for (i = 0; i < sn; i++) - D(i) = fix_mul(D(i), 10078); /* 10076 */ - for (i = 0; i < dn; i++) - S(i) = fix_mul(S(i), 13318); /* 13320 */ - for (i = 0; i < sn; i++) - D(i) -= fix_mul(SS_(i) + SS_(i + 1), 3633); - for (i = 0; i < dn; i++) - S(i) -= fix_mul(DD_(i) + DD_(i - 1), 7233); - for (i = 0; i < sn; i++) - D(i) += fix_mul(SS_(i) + SS_(i + 1), 434); - for (i = 0; i < dn; i++) - S(i) += fix_mul(DD_(i) + DD_(i - 1), 12993); - } - } else { + int i; + if (!cas) { if ((dn > 0) || (sn > 1)) { /* NEW : CASE ONE ELEMENT */ for (i = 0; i < sn; i++) - S(i) = fix_mul(S(i), 10078); /* 10076 */ + S(i) = fix_mul(S(i), 10078); /* 10076 */ for (i = 0; i < dn; i++) - D(i) = fix_mul(D(i), 13318); /* 13320 */ + D(i) = fix_mul(D(i), 13318); /* 13320 */ for (i = 0; i < sn; i++) - S(i) -= fix_mul(D_(i - 1) + D_(i), 3633); + S(i) -= fix_mul(D_(i - 1) + D_(i), 3633); for (i = 0; i < dn; i++) - D(i) -= fix_mul(S_(i) + S_(i + 1), 7233); + D(i) -= fix_mul(S_(i) + S_(i + 1), 7233); for (i = 0; i < sn; i++) - S(i) += fix_mul(D_(i - 1) + D_(i), 434); + S(i) += fix_mul(D_(i - 1) + D_(i), 434); for (i = 0; i < dn; i++) - D(i) += fix_mul(S_(i) + S_(i + 1), 12993); + D(i) += fix_mul(S_(i) + S_(i + 1), 12993); + } + } else { + if ((sn > 0) || (dn > 1)) { /* NEW : CASE ONE ELEMENT */ + for (i = 0; i < sn; i++) + D(i) = fix_mul(D(i), 10078); /* 10076 */ + for (i = 0; i < dn; i++) + S(i) = fix_mul(S(i), 13318); /* 13320 */ + for (i = 0; i < sn; i++) + D(i) -= fix_mul(SS_(i) + SS_(i + 1), 3633); + for (i = 0; i < dn; i++) + S(i) -= fix_mul(DD_(i) + DD_(i - 1), 7233); + for (i = 0; i < sn; i++) + D(i) += fix_mul(SS_(i) + SS_(i + 1), 434); + for (i = 0; i < dn; i++) + S(i) += fix_mul(DD_(i) + DD_(i - 1), 12993); } } } @@ -376,67 +412,118 @@ void dwt_decode_1_real(int *a, int n, int x, int res, int cas) /* Forward 9-7 wavelet transform in 2-D. */ /* */ -void dwt_encode_real(int *a, int w, int h, tcd_tilecomp_t * tilec, int l) +void dwt_encode_real(tcd_tilecomp_t * tilec) { - int i, j; - int rw; /* width of the resolution level computed */ - int rh; /* heigth of the resolution level computed */ - int rw1; /* width of the resolution level once lower than computed one */ - int rh1; /* height of the resolution level once lower than computed one */ + int i, j, k; + int* a; + int* aj; + int* bj; + int w, l; + + w = tilec->x1-tilec->x0; + l = tilec->numresolutions-1; + a = tilec->data; for (i = 0; i < l; i++) { - int cas_col = 0; /* 0 = non inversion on horizontal filtering 1 = inversion between low-pass and high-pass filtering */ - int cas_row = 0; /* 0 = non inversion on vertical filtering 1 = inversion between low-pass and high-pass filtering */ + int rw; /* width of the resolution level computed */ + int rh; /* heigth of the resolution level computed */ + int rw1; /* width of the resolution level once lower than computed one */ + int rh1; /* height of the resolution level once lower than computed one */ + int cas_col; /* 0 = non inversion on horizontal filtering 1 = inversion between low-pass and high-pass filtering */ + int cas_row; /* 0 = non inversion on vertical filtering 1 = inversion between low-pass and high-pass filtering */ + int dn, sn; + rw = tilec->resolutions[l - i].x1 - tilec->resolutions[l - i].x0; rh = tilec->resolutions[l - i].y1 - tilec->resolutions[l - i].y0; - rw1 = - tilec->resolutions[l - i - 1].x1 - tilec->resolutions[l - i - 1].x0; - rh1 = - tilec->resolutions[l - i - 1].y1 - tilec->resolutions[l - i - 1].y0; + rw1= tilec->resolutions[l - i - 1].x1 - tilec->resolutions[l - i - 1].x0; + rh1= tilec->resolutions[l - i - 1].y1 - tilec->resolutions[l - i - 1].y0; cas_row = tilec->resolutions[l - i].x0 % 2; cas_col = tilec->resolutions[l - i].y0 % 2; - for (j = 0; j < rw; j++) - dwt_encode_1_real(a + j, rh, w, rh1, cas_col); - for (j = 0; j < rh; j++) - dwt_encode_1_real(a + j * w, rw, 1, rw1, cas_row); + sn = rh1; + dn = rh - rh1; + bj=(int*)malloc(rh*sizeof(int)); + for (j = 0; j < rw; j++) { + aj = a + j; + for (k = 0; k < rh; k++) bj[k] = aj[k*w]; + dwt_encode_1_real(bj, dn, sn, cas_col); + dwt_deinterleave_v(bj, aj, dn, sn, w, cas_col); + } + free(bj); + + sn = rw1; + dn = rw - rw1; + bj=(int*)malloc(rw*sizeof(int)); + for (j = 0; j < rh; j++) { + aj = a + j * w; + for (k = 0; k < rw; k++) bj[k] = aj[k]; + dwt_encode_1_real(bj, dn, sn, cas_row); + dwt_deinterleave_h(bj, aj, dn, sn, cas_row); + } + free(bj); } } + /* */ /* Inverse 9-7 wavelet transform in 2-D. */ /* */ -void dwt_decode_real(int *a, int w, int h, tcd_tilecomp_t * tilec, int l, - int stop) +void dwt_decode_real(tcd_tilecomp_t * tilec, int stop) { - int i, j; - int rw; /* width of the resolution level computed */ - int rh; /* heigth of the resolution level computed */ - int rw1; /* width of the resolution level once lower than computed one */ - int rh1; /* height of the resolution level once lower than computed one */ - for (i = l - 1; i >= stop; i--) { - int cas_col = 0; /* 0 = non inversion on horizontal filtering 1 = inversion between low-pass and high-pass filtering */ - int cas_row = 0; /* 0 = non inversion on vertical filtering 1 = inversion between low-pass and high-pass filtering */ + int i, j, k; + int* a; + int* aj; + int* bj; + int w, l; + + w = tilec->x1-tilec->x0; + l = tilec->numresolutions-1; + a = tilec->data; + + for (i = l-1; i >= stop; i--) { + int rw; /* width of the resolution level computed */ + int rh; /* heigth of the resolution level computed */ + int rw1; /* width of the resolution level once lower than computed one */ + int rh1; /* height of the resolution level once lower than computed one */ + int cas_col; /* 0 = non inversion on horizontal filtering 1 = inversion between low-pass and high-pass filtering */ + int cas_row; /* 0 = non inversion on vertical filtering 1 = inversion between low-pass and high-pass filtering */ + int dn, sn; rw = tilec->resolutions[l - i].x1 - tilec->resolutions[l - i].x0; rh = tilec->resolutions[l - i].y1 - tilec->resolutions[l - i].y0; - rw1 = - tilec->resolutions[l - i - 1].x1 - tilec->resolutions[l - i - 1].x0; - rh1 = - tilec->resolutions[l - i - 1].y1 - tilec->resolutions[l - i - 1].y0; + rw1= tilec->resolutions[l - i - 1].x1 - tilec->resolutions[l - i - 1].x0; + rh1= tilec->resolutions[l - i - 1].y1 - tilec->resolutions[l - i - 1].y0; - cas_row = tilec->resolutions[l - i].x0 % 2; - cas_col = tilec->resolutions[l - i].y0 % 2; + cas_col = tilec->resolutions[l - i].x0 % 2; /* 0 = non inversion on horizontal filtering 1 = inversion between low-pass and high-pass filtering */ + cas_row = tilec->resolutions[l - i].y0 % 2; /* 0 = non inversion on vertical filtering 1 = inversion between low-pass and high-pass filtering */ - for (j = 0; j < rh; j++) - dwt_decode_1_real(a + j * w, rw, 1, rw1, cas_row); - for (j = 0; j < rw; j++) - dwt_decode_1_real(a + j, rh, w, rh1, cas_col); + sn = rw1; + dn = rw-rw1; + bj = (int*)malloc(rw * sizeof(int)); + for (j = 0; j < rh; j++) { + aj = a+j*w; + dwt_interleave_h(aj, bj, dn, sn, cas_col); + dwt_decode_1_real(bj, dn, sn, cas_col); + for (k = 0; k < rw; k++) aj[k] = bj[k]; + } + free(bj); + + sn = rh1; + dn = rh-rh1; + bj = (int*)malloc(rh * sizeof(int)); + for (j=0; j */ /* Get gain of 9-7 wavelet transform. */ /* */ diff --git a/libopenjpeg/dwt.h b/libopenjpeg/dwt.h index ec9860b6..a9c30ca7 100644 --- a/libopenjpeg/dwt.h +++ b/libopenjpeg/dwt.h @@ -33,26 +33,15 @@ /* * Apply a reversible DWT transform to a component of an image - * a: samples of the component - * w: width of the component - * h: height of the component * tilec : tile component information (present tile) - * l: number of decomposition levels in the DWT */ /* void dwt_encode(int* a, int w, int h, int l); */ -void dwt_encode(int *a, int w, int h, tcd_tilecomp_t * tilec, int l); +void dwt_encode(tcd_tilecomp_t * tilec); /* * Apply a reversible inverse DWT transform to a component of an image - * a: samples of the component - * w: width of the component - * h: height of the component * tilec : tile component information (present tile) - * l: number of decomposition levels in the DWT - * row_tilec : tile component information (previous tile on the same row) - * col_tilec : tile component information (previous tile on the same column) */ -void dwt_decode(int *a, int w, int h, tcd_tilecomp_t * tilec, int l, - int stop); +void dwt_decode(tcd_tilecomp_t * tilec, int stop); /* * Get the gain of a subband for the reversible DWT @@ -69,22 +58,13 @@ double dwt_getnorm(int level, int orient); /* * Apply an irreversible DWT transform to a component of an image - * a: samples of the component - * w: width of the component - * h: height of the component - * l: number of decomposition levels in the DWT */ -void dwt_encode_real(int *a, int w, int h, tcd_tilecomp_t * tilec, int l); +void dwt_encode_real(tcd_tilecomp_t * tilec); /* * Apply an irreversible inverse DWT transform to a component of an image - * a: samples of the component - * w: width of the component - * h: height of the component - * l: number of decomposition levels in the DWT */ -void dwt_decode_real(int *a, int w, int h, tcd_tilecomp_t * tilec, int l, - int stop); +void dwt_decode_real(tcd_tilecomp_t * tilec, int stop); /* * Get the gain of a subband for the irreversible DWT * orient: number that identifies the subband (0->LL, 1->HL, 2->LH, 3->HH) diff --git a/libopenjpeg/tcd.c b/libopenjpeg/tcd.c index ec1f6968..3e433209 100644 --- a/libopenjpeg/tcd.c +++ b/libopenjpeg/tcd.c @@ -1161,7 +1161,7 @@ tcd_encode_tile_pxm(int tileno, unsigned char *dest, int len, tcd_tile_t *tile; j2k_tcp_t *tcp = &tcd_cp->tcps[0]; j2k_tccp_t *tccp = &tcp->tccps[0]; - + tcd_tileno = tileno; tcd_tile = tcd_image.tiles; tcd_tcp = &tcd_cp->tcps[tileno]; @@ -1169,19 +1169,19 @@ tcd_encode_tile_pxm(int tileno, unsigned char *dest, int len, /* INDEX >> "Precinct_nb_X et Precinct_nb_Y" */ if (info_IM->index_on) { tcd_tilecomp_t *tilec_idx = &tile->comps[0]; //Based on Component 0 - + for (i = 0; i < tilec_idx->numresolutions; i++) { - + tcd_resolution_t *res_idx = &tilec_idx->resolutions[i]; - + info_IM->tile[tileno].pw[i] = res_idx->pw; info_IM->tile[tileno].ph[i] = res_idx->ph; npck+=res_idx->pw * res_idx->ph; - + info_IM->tile[tileno].pdx[i] = tccp->prcw[i]; info_IM->tile[tileno].pdy[i] = tccp->prch[i]; - + } info_IM->tile[tileno].packet = (info_packet *) calloc(info_IM->Comp * info_IM->Layer * npck, sizeof(info_packet)); } @@ -1261,18 +1261,16 @@ tcd_encode_tile_pxm(int tileno, unsigned char *dest, int len, } /*----------------DWT---------------------*/ - /* time3=clock(); */ - for (compno = 0; compno < tile->numcomps; compno++) { - tcd_tilecomp_t *tilec = &tile->comps[compno]; - if (tcd_tcp->tccps[compno].qmfbid == 1) { - dwt_encode(tilec->data, tilec->x1 - tilec->x0, - tilec->y1 - tilec->y0, tilec, tilec->numresolutions - 1); - } else if (tcd_tcp->tccps[compno].qmfbid == 0) { - dwt_encode_real(tilec->data, tilec->x1 - tilec->x0, - tilec->y1 - tilec->y0, tilec, - tilec->numresolutions - 1); - } +// mod Ive +for (compno = 0; compno < tile->numcomps; compno++) { + tcd_tilecomp_t *tilec = &tile->comps[compno]; + if (tcd_tcp->tccps[compno].qmfbid == 1) { + dwt_encode(tilec); + } else if (tcd_tcp->tccps[compno].qmfbid == 0) { + dwt_encode_real(tilec); } +} +// /mod Ive /*------------------TIER1-----------------*/ t1_init_luts(); @@ -1317,7 +1315,7 @@ tcd_encode_tile_pgx(int tileno, unsigned char *dest, int len, tcd_tile_t *tile; j2k_tcp_t *tcp = &tcd_cp->tcps[0]; j2k_tccp_t *tccp = &tcp->tccps[0]; - + tcd_tileno = tileno; tcd_tile = tcd_image.tiles; tcd_tcp = &tcd_cp->tcps[tileno]; @@ -1415,17 +1413,16 @@ tcd_encode_tile_pgx(int tileno, unsigned char *dest, int len, /*----------------DWT---------------------*/ - for (compno = 0; compno < tile->numcomps; compno++) { - tcd_tilecomp_t *tilec = &tile->comps[compno]; - if (tcd_tcp->tccps[compno].qmfbid == 1) { - dwt_encode(tilec->data, tilec->x1 - tilec->x0, - tilec->y1 - tilec->y0, tilec, tilec->numresolutions - 1); - } else if (tcd_tcp->tccps[compno].qmfbid == 0) { - dwt_encode_real(tilec->data, tilec->x1 - tilec->x0, - tilec->y1 - tilec->y0, tilec, - tilec->numresolutions - 1); - } +// mod Ive +for (compno = 0; compno < tile->numcomps; compno++) { + tcd_tilecomp_t *tilec = &tile->comps[compno]; + if (tcd_tcp->tccps[compno].qmfbid == 1) { + dwt_encode(tilec); + } else if (tcd_tcp->tccps[compno].qmfbid == 0) { + dwt_encode_real(tilec); } +} +// /mod Ive /*------------------TIER1-----------------*/ @@ -1508,20 +1505,18 @@ int tcd_decode_tile(unsigned char *src, int len, int tileno) } + // mod Ive if (tcd_tcp->tccps[compno].qmfbid == 1) { - dwt_decode(tilec->data, tilec->x1 - tilec->x0, - tilec->y1 - tilec->y0, tilec, - tilec->numresolutions - 1, - tilec->numresolutions - 1 - - tcd_img->comps[compno].resno_decoded); + dwt_decode(tilec, + tilec->numresolutions - 1 - + tcd_img->comps[compno].resno_decoded); } else { - dwt_decode_real(tilec->data, tilec->x1 - tilec->x0, - tilec->y1 - tilec->y0, tilec, - tilec->numresolutions - 1, + dwt_decode_real(tilec, tilec->numresolutions - 1 - tcd_img->comps[compno].resno_decoded); } - + // /mod Ive + if (tile->comps[compno].numresolutions > 0) tcd_img->comps[compno].factor = tile->comps[compno].numresolutions -