walkero
/
openjpeg
1
0
Fork 0
Code Issues Pull Requests Projects Releases 1 Wiki Activity

Irreversible compression/decompression DWT: use 1/K constant as per standard 

The previous constant opj_c13318 was mysteriously equal to 2/K , and in
the DWT, we had to divide K and opj_c13318 by 2... The issue was that the
band->stepsize computation in tcd.c didn't take into account the log2gain of
the band.

The effect of this change is expected to be mostly equivalent to the previous
situation, except some difference in rounding. But it leads to a dramatic
reduction of the mean square error and peak error in the irreversible encoding
of issue141.tif !
This commit is contained in:
Even Rouault 2020-05-20 18:00:45 +02:00
parent f38c069547
commit 3cd1305596
No known key found for this signature in database
GPG Key ID: 33EBBFC47B3DD87D
5 changed files with 33 additions and 67 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

44
src/lib/openjp2/dwt.c
View File

@ -103,13 +103,13 @@ typedef struct v4dwt_local {
} opj_v4dwt_t ; } opj_v4dwt_t ;
/* From table F.4 from the standard */ /* From table F.4 from the standard */
static const OPJ_FLOAT32 opj_dwt_alpha = -1.586134342f; /* 12994 */ static const OPJ_FLOAT32 opj_dwt_alpha = -1.586134342f;
static const OPJ_FLOAT32 opj_dwt_beta = -0.052980118f; /* 434 */ static const OPJ_FLOAT32 opj_dwt_beta = -0.052980118f;
static const OPJ_FLOAT32 opj_dwt_gamma = 0.882911075f; /* -7233 */ static const OPJ_FLOAT32 opj_dwt_gamma = 0.882911075f;
static const OPJ_FLOAT32 opj_dwt_delta = 0.443506852f; /* -3633 */ static const OPJ_FLOAT32 opj_dwt_delta = 0.443506852f;
static const OPJ_FLOAT32 opj_K = 1.230174105f; /* 10078 */ static const OPJ_FLOAT32 opj_K = 1.230174105f;
static const OPJ_FLOAT32 opj_c13318 = 1.625732422f; static const OPJ_FLOAT32 opj_invK = (OPJ_FLOAT32)(1.0 / 1.230174105);
/*@}*/ /*@}*/
@ -1108,9 +1108,9 @@ static void opj_dwt_encode_1_real(void *aIn, OPJ_INT32 dn, OPJ_INT32 sn,
(OPJ_UINT32)opj_int_min(sn, dn - a), (OPJ_UINT32)opj_int_min(sn, dn - a),
opj_dwt_delta); opj_dwt_delta);
opj_dwt_encode_step1(w + b, 0, (OPJ_UINT32)dn, opj_dwt_encode_step1(w + b, 0, (OPJ_UINT32)dn,
opj_K / 2); opj_K);
opj_dwt_encode_step1(w + a, 0, (OPJ_UINT32)sn, opj_dwt_encode_step1(w + a, 0, (OPJ_UINT32)sn,
opj_c13318 / 2); opj_invK);
} }
static void opj_dwt_encode_stepsize(OPJ_INT32 stepsize, OPJ_INT32 numbps, static void opj_dwt_encode_stepsize(OPJ_INT32 stepsize, OPJ_INT32 numbps,
@ -1399,21 +1399,6 @@ OPJ_BOOL opj_dwt_decode(opj_tcd_t *p_tcd, opj_tcd_tilecomp_t* tilec,
} }
} }
/* <summary> */
/* Get gain of 5-3 wavelet transform. */
/* </summary> */
OPJ_UINT32 opj_dwt_getgain(OPJ_UINT32 orient)
{
if (orient == 0) {
return 0;
}
if (orient == 1 || orient == 2) {
return 1;
}
return 2;
}
/* <summary> */ /* <summary> */
/* Get norm of 5-3 wavelet. */ /* Get norm of 5-3 wavelet. */
/* </summary> */ /* </summary> */
@ -1440,15 +1425,6 @@ OPJ_BOOL opj_dwt_encode_real(opj_tcd_t *p_tcd,
opj_dwt_encode_1_real); opj_dwt_encode_1_real);
} }
/* <summary> */
/* Get gain of 9-7 wavelet transform. */
/* </summary> */
OPJ_UINT32 opj_dwt_getgain_real(OPJ_UINT32 orient)
{
(void)orient;
return 0;
}
/* <summary> */ /* <summary> */
/* Get norm of 9-7 wavelet. */ /* Get norm of 9-7 wavelet. */
/* </summary> */ /* </summary> */
@ -2649,7 +2625,7 @@ static void opj_v4dwt_decode(opj_v4dwt_t* OPJ_RESTRICT dwt)
opj_v4dwt_decode_step1_sse(dwt->wavelet + a, dwt->win_l_x0, dwt->win_l_x1, opj_v4dwt_decode_step1_sse(dwt->wavelet + a, dwt->win_l_x0, dwt->win_l_x1,
_mm_set1_ps(opj_K)); _mm_set1_ps(opj_K));
opj_v4dwt_decode_step1_sse(dwt->wavelet + b, dwt->win_h_x0, dwt->win_h_x1, opj_v4dwt_decode_step1_sse(dwt->wavelet + b, dwt->win_h_x0, dwt->win_h_x1,
_mm_set1_ps(opj_c13318)); _mm_set1_ps(opj_invK));
opj_v4dwt_decode_step2_sse(dwt->wavelet + b, dwt->wavelet + a + 1, opj_v4dwt_decode_step2_sse(dwt->wavelet + b, dwt->wavelet + a + 1,
dwt->win_l_x0, dwt->win_l_x1, dwt->win_l_x0, dwt->win_l_x1,
(OPJ_UINT32)opj_int_min(dwt->sn, dwt->dn - a), (OPJ_UINT32)opj_int_min(dwt->sn, dwt->dn - a),
@ -2670,7 +2646,7 @@ static void opj_v4dwt_decode(opj_v4dwt_t* OPJ_RESTRICT dwt)
opj_v4dwt_decode_step1(dwt->wavelet + a, dwt->win_l_x0, dwt->win_l_x1, opj_v4dwt_decode_step1(dwt->wavelet + a, dwt->win_l_x0, dwt->win_l_x1,
opj_K); opj_K);
opj_v4dwt_decode_step1(dwt->wavelet + b, dwt->win_h_x0, dwt->win_h_x1, opj_v4dwt_decode_step1(dwt->wavelet + b, dwt->win_h_x0, dwt->win_h_x1,
opj_c13318); opj_invK);
opj_v4dwt_decode_step2(dwt->wavelet + b, dwt->wavelet + a + 1, opj_v4dwt_decode_step2(dwt->wavelet + b, dwt->wavelet + a + 1,
dwt->win_l_x0, dwt->win_l_x1, dwt->win_l_x0, dwt->win_l_x1,
(OPJ_UINT32)opj_int_min(dwt->sn, dwt->dn - a), (OPJ_UINT32)opj_int_min(dwt->sn, dwt->dn - a),

12
src/lib/openjp2/dwt.h
View File

@ -73,12 +73,6 @@ OPJ_BOOL opj_dwt_decode(opj_tcd_t *p_tcd,
opj_tcd_tilecomp_t* tilec, opj_tcd_tilecomp_t* tilec,
OPJ_UINT32 numres); OPJ_UINT32 numres);
/**
Get the gain of a subband for the reversible 5-3 DWT.
@param orient Number that identifies the subband (0->LL, 1->HL, 2->LH, 3->HH)
@return Returns 0 if orient = 0, returns 1 if orient = 1 or 2, returns 2 otherwise
*/
OPJ_UINT32 opj_dwt_getgain(OPJ_UINT32 orient) ;
/** /**
Get the norm of a wavelet function of a subband at a specified level for the reversible 5-3 DWT. Get the norm of a wavelet function of a subband at a specified level for the reversible 5-3 DWT.
@param level Level of the wavelet function @param level Level of the wavelet function
@ -105,12 +99,6 @@ OPJ_BOOL opj_dwt_decode_real(opj_tcd_t *p_tcd,
opj_tcd_tilecomp_t* OPJ_RESTRICT tilec, opj_tcd_tilecomp_t* OPJ_RESTRICT tilec,
OPJ_UINT32 numres); OPJ_UINT32 numres);
/**
Get the gain of a subband for the irreversible 9-7 DWT.
@param orient Number that identifies the subband (0->LL, 1->HL, 2->LH, 3->HH)
@return Returns the gain of the 9-7 wavelet transform
*/
OPJ_UINT32 opj_dwt_getgain_real(OPJ_UINT32 orient);
/** /**
Get the norm of a wavelet function of a subband at a specified level for the irreversible 9-7 DWT Get the norm of a wavelet function of a subband at a specified level for the irreversible 9-7 DWT
@param level Level of the wavelet function @param level Level of the wavelet function

4
src/lib/openjp2/t1.c
View File

@ -1426,7 +1426,11 @@ static OPJ_FLOAT64 opj_t1_getwmsedec(
if (qmfbid == 1) { if (qmfbid == 1) {
w2 = opj_dwt_getnorm(level, orient); w2 = opj_dwt_getnorm(level, orient);
} else { /* if (qmfbid == 0) */ } else { /* if (qmfbid == 0) */
const OPJ_INT32 log2_gain = (orient == 0) ? 0 :
(orient == 3) ? 2 : 1;
w2 = opj_dwt_getnorm_real(level, orient); w2 = opj_dwt_getnorm_real(level, orient);
/* Not sure this is right. But preserves past behaviour */
stepsize /= (1 << log2_gain);
} }
wmsedec = w1 * w2 * stepsize * (1 << bpno); wmsedec = w1 * w2 * stepsize * (1 << bpno);

28
src/lib/openjp2/tcd.c
View File

@ -724,7 +724,6 @@ static INLINE OPJ_BOOL opj_tcd_init_tile(opj_tcd_t *p_tcd, OPJ_UINT32 p_tile_no,
OPJ_BOOL isEncoder, OPJ_SIZE_T sizeof_block, OPJ_BOOL isEncoder, OPJ_SIZE_T sizeof_block,
opj_event_mgr_t* manager) opj_event_mgr_t* manager)
{ {
OPJ_UINT32(*l_gain_ptr)(OPJ_UINT32) = 00;
OPJ_UINT32 compno, resno, bandno, precno, cblkno; OPJ_UINT32 compno, resno, bandno, precno, cblkno;
opj_tcp_t * l_tcp = 00; opj_tcp_t * l_tcp = 00;
opj_cp_t * l_cp = 00; opj_cp_t * l_cp = 00;
@ -740,7 +739,6 @@ static INLINE OPJ_BOOL opj_tcd_init_tile(opj_tcd_t *p_tcd, OPJ_UINT32 p_tile_no,
OPJ_UINT32 p, q; OPJ_UINT32 p, q;
OPJ_UINT32 l_level_no; OPJ_UINT32 l_level_no;
OPJ_UINT32 l_pdx, l_pdy; OPJ_UINT32 l_pdx, l_pdy;
OPJ_UINT32 l_gain;
OPJ_INT32 l_x0b, l_y0b; OPJ_INT32 l_x0b, l_y0b;
OPJ_UINT32 l_tx0, l_ty0; OPJ_UINT32 l_tx0, l_ty0;
/* extent of precincts , top left, bottom right**/ /* extent of precincts , top left, bottom right**/
@ -879,11 +877,6 @@ static INLINE OPJ_BOOL opj_tcd_init_tile(opj_tcd_t *p_tcd, OPJ_UINT32 p_tile_no,
l_level_no = l_tilec->numresolutions; l_level_no = l_tilec->numresolutions;
l_res = l_tilec->resolutions; l_res = l_tilec->resolutions;
l_step_size = l_tccp->stepsizes; l_step_size = l_tccp->stepsizes;
if (l_tccp->qmfbid == 0) {
l_gain_ptr = &opj_dwt_getgain_real;
} else {
l_gain_ptr = &opj_dwt_getgain;
}
/*fprintf(stderr, "\tlevel_no=%d\n",l_level_no);*/ /*fprintf(stderr, "\tlevel_no=%d\n",l_level_no);*/
for (resno = 0; resno < l_tilec->numresolutions; ++resno) { for (resno = 0; resno < l_tilec->numresolutions; ++resno) {
@ -970,7 +963,6 @@ static INLINE OPJ_BOOL opj_tcd_init_tile(opj_tcd_t *p_tcd, OPJ_UINT32 p_tile_no,
l_band = l_res->bands; l_band = l_res->bands;
for (bandno = 0; bandno < l_res->numbands; ++bandno, ++l_band, ++l_step_size) { for (bandno = 0; bandno < l_res->numbands; ++bandno, ++l_band, ++l_step_size) {
OPJ_INT32 numbps;
/*fprintf(stderr, "\t\t\tband_no=%d/%d\n", bandno, l_res->numbands );*/ /*fprintf(stderr, "\t\t\tband_no=%d/%d\n", bandno, l_res->numbands );*/
if (resno == 0) { if (resno == 0) {
@ -1006,14 +998,20 @@ static INLINE OPJ_BOOL opj_tcd_init_tile(opj_tcd_t *p_tcd, OPJ_UINT32 p_tile_no,
} }
} }
/** avoid an if with storing function pointer */ {
l_gain = (*l_gain_ptr)(l_band->bandno); /* Table E-1 - Sub-band gains */
numbps = (OPJ_INT32)(l_image_comp->prec + l_gain); const OPJ_INT32 log2_gain = (l_band->bandno == 0) ? 0 :
(l_band->bandno == 3) ? 2 : 1;
/* Nominal dynamic range. Equation E-4 */
const OPJ_INT32 Rb = (OPJ_INT32)l_image_comp->prec + log2_gain;
/* Delta_b value of Equation E-3 in "E.1 Inverse quantization
* procedure" of the standard */
l_band->stepsize = (OPJ_FLOAT32)(((1.0 + l_step_size->mant / 2048.0) * pow(2.0,
(OPJ_INT32)(Rb - l_step_size->expn))));
}
/* Delta_b value of Equation E-3 in "E.1 Inverse quantization
* procedure" of the standard */
l_band->stepsize = (OPJ_FLOAT32)(((1.0 + l_step_size->mant / 2048.0) * pow(2.0,
(OPJ_INT32)(numbps - l_step_size->expn))));
/* Mb value of Equation E-2 in "E.1 Inverse quantization /* Mb value of Equation E-2 in "E.1 Inverse quantization
* procedure" of the standard */ * procedure" of the standard */
l_band->numbps = l_step_size->expn + (OPJ_INT32)l_tccp->numgbits - l_band->numbps = l_step_size->expn + (OPJ_INT32)l_tccp->numgbits -

12
tests/nonregression/test_suite.ctest.in
View File

@ -32,16 +32,16 @@ opj_compress -i @INPUT_NR_PATH@/random-issue-0005.tif -o @TEMP_PATH@/random-issu
# related to issue 62 # related to issue 62
opj_compress -i @INPUT_NR_PATH@/tmp-issue-0062.raw -o @TEMP_PATH@/tmp-issue-0062-u.raw.j2k -F 512,512,1,16,u opj_compress -i @INPUT_NR_PATH@/tmp-issue-0062.raw -o @TEMP_PATH@/tmp-issue-0062-u.raw.j2k -F 512,512,1,16,u
opj_compress -i @INPUT_NR_PATH@/tmp-issue-0062.raw -o @TEMP_PATH@/tmp-issue-0062-s.raw.j2k -F 512,512,1,16,s opj_compress -i @INPUT_NR_PATH@/tmp-issue-0062.raw -o @TEMP_PATH@/tmp-issue-0062-s.raw.j2k -F 512,512,1,16,s
opj_compress lossy-check { -n 3 -i prec -m 175:100:212 -p 78:63:91 } -i @INPUT_NR_PATH@/X_4_2K_24_185_CBR_WB_000.tif -o @TEMP_PATH@/X_4_2K_24_185_CBR_WB_000_C2K_24.j2k -cinema2K 24 opj_compress lossy-check { -n 3 -i prec -m 175:100:212 -p 79:64:92 } -i @INPUT_NR_PATH@/X_4_2K_24_185_CBR_WB_000.tif -o @TEMP_PATH@/X_4_2K_24_185_CBR_WB_000_C2K_24.j2k -cinema2K 24
opj_compress lossy-check { -n 3 -i prec -m 298:168:363 -p 121:73:164 } -i @INPUT_NR_PATH@/X_5_2K_24_235_CBR_STEM24_000.tif -o @TEMP_PATH@/X_5_2K_24_235_CBR_STEM24_000_C2K_24.j2k -cinema2K 24 opj_compress lossy-check { -n 3 -i prec -m 298:168:363 -p 122:73:164 } -i @INPUT_NR_PATH@/X_5_2K_24_235_CBR_STEM24_000.tif -o @TEMP_PATH@/X_5_2K_24_235_CBR_STEM24_000_C2K_24.j2k -cinema2K 24
opj_compress lossy-check { -n 3 -i prec -m 76:54:140 -p 55:49:74 } -i @INPUT_NR_PATH@/X_6_2K_24_FULL_CBR_CIRCLE_000.tif -o @TEMP_PATH@/X_6_2K_24_FULL_CBR_CIRCLE_000_C2K_24.j2k -cinema2K 24 opj_compress lossy-check { -n 3 -i prec -m 76:54:140 -p 56:49:74 } -i @INPUT_NR_PATH@/X_6_2K_24_FULL_CBR_CIRCLE_000.tif -o @TEMP_PATH@/X_6_2K_24_FULL_CBR_CIRCLE_000_C2K_24.j2k -cinema2K 24
opj_compress lossy-check { -n 3 -i prec -m 384:385:842 -p 134:146:200 } -i @INPUT_NR_PATH@/X_4_2K_24_185_CBR_WB_000.tif -o @TEMP_PATH@/X_4_2K_24_185_CBR_WB_000_C2K_48.j2k -cinema2K 48 opj_compress lossy-check { -n 3 -i prec -m 384:385:842 -p 135:144:202 } -i @INPUT_NR_PATH@/X_4_2K_24_185_CBR_WB_000.tif -o @TEMP_PATH@/X_4_2K_24_185_CBR_WB_000_C2K_48.j2k -cinema2K 48
opj_compress lossy-check { -n 3 -i prec -m 933:827:2206 -p 201:184:314 } -i @INPUT_NR_PATH@/X_5_2K_24_235_CBR_STEM24_000.tif -o @TEMP_PATH@/X_5_2K_24_235_CBR_STEM24_000_C2K_48.j2k -cinema2K 48 opj_compress lossy-check { -n 3 -i prec -m 933:827:2206 -p 201:184:314 } -i @INPUT_NR_PATH@/X_5_2K_24_235_CBR_STEM24_000.tif -o @TEMP_PATH@/X_5_2K_24_235_CBR_STEM24_000_C2K_48.j2k -cinema2K 48
opj_compress lossy-check { -n 3 -i prec -m 194:173:531 -p 94:79:154 } -i @INPUT_NR_PATH@/X_6_2K_24_FULL_CBR_CIRCLE_000.tif -o @TEMP_PATH@/X_6_2K_24_FULL_CBR_CIRCLE_000_C2K_48.j2k -cinema2K 48 opj_compress lossy-check { -n 3 -i prec -m 194:173:531 -p 94:79:154 } -i @INPUT_NR_PATH@/X_6_2K_24_FULL_CBR_CIRCLE_000.tif -o @TEMP_PATH@/X_6_2K_24_FULL_CBR_CIRCLE_000_C2K_48.j2k -cinema2K 48
opj_compress lossy-check { -n 3 -i prec -m 6:4:7 -p 141:141:193 } -i @INPUT_NR_PATH@/ElephantDream_4K.tif -o @TEMP_PATH@/ElephantDream_4K_C4K.j2k -cinema4K opj_compress lossy-check { -n 3 -i prec -m 6:4:7 -p 141:141:191 } -i @INPUT_NR_PATH@/ElephantDream_4K.tif -o @TEMP_PATH@/ElephantDream_4K_C4K.j2k -cinema4K
# issue 141 # issue 141
opj_compress -i @INPUT_NR_PATH@/issue141.rawl -o @TEMP_PATH@/issue141.rawl.j2k -F 2048,32,1,16,u opj_compress -i @INPUT_NR_PATH@/issue141.rawl -o @TEMP_PATH@/issue141.rawl.j2k -F 2048,32,1,16,u
opj_compress lossy-check { -n 1 -m 61 -p 11 } -i @INPUT_NR_PATH@/issue141.tif -o @TEMP_PATH@/issue141-I.rawl.j2k -I opj_compress lossy-check { -n 1 -m 0.1 -p 2 } -i @INPUT_NR_PATH@/issue141.tif -o @TEMP_PATH@/issue141-I.rawl.j2k -I
# issue 46: # issue 46:
opj_compress -i @INPUT_NR_PATH@/Bretagne2.ppm -o @TEMP_PATH@/Bretagne2_5.j2k -c [64,64] opj_compress -i @INPUT_NR_PATH@/Bretagne2.ppm -o @TEMP_PATH@/Bretagne2_5.j2k -c [64,64]
# issue 316 # issue 316