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 !

diff --git a/src/lib/openjp2/dwt.c b/src/lib/openjp2/dwt.c
index a825f0136082c8b94194228aa2ee80ca847950dd..de8fdf4e9ea28b3dd5a79ba24dd8555df14219cc 100644 (file)
--- a/src/lib/openjp2/dwt.c
+++ b/src/lib/openjp2/dwt.c
@@ -103,13 +103,13 @@ typedef struct v4dwt_local {
 } opj_v4dwt_t ;
 
 /* From table F.4 from the standard */
-static const OPJ_FLOAT32 opj_dwt_alpha =  -1.586134342f; /*  12994 */
-static const OPJ_FLOAT32 opj_dwt_beta  =  -0.052980118f; /*    434 */
-static const OPJ_FLOAT32 opj_dwt_gamma = 0.882911075f; /*  -7233 */
-static const OPJ_FLOAT32 opj_dwt_delta = 0.443506852f; /*  -3633 */
+static const OPJ_FLOAT32 opj_dwt_alpha =  -1.586134342f;
+static const OPJ_FLOAT32 opj_dwt_beta  =  -0.052980118f;
+static const OPJ_FLOAT32 opj_dwt_gamma = 0.882911075f;
+static const OPJ_FLOAT32 opj_dwt_delta = 0.443506852f;
 
-static const OPJ_FLOAT32 opj_K      = 1.230174105f; /*  10078 */
-static const OPJ_FLOAT32 opj_c13318 = 1.625732422f;
+static const OPJ_FLOAT32 opj_K      = 1.230174105f;
+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_dwt_delta);
     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_c13318 / 2);
+                         opj_invK);
 }
 
 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>                */
 /* Get norm of 5-3 wavelet. */
 /* </summary>               */
@@ -1440,15 +1425,6 @@ OPJ_BOOL opj_dwt_encode_real(opj_tcd_t *p_tcd,
                                     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>                */
 /* Get norm of 9-7 wavelet. */
 /* </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,
                                _mm_set1_ps(opj_K));
     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,
                                dwt->win_l_x0, dwt->win_l_x1,
                                (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_K);
     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,
                            dwt->win_l_x0, dwt->win_l_x1,
                            (OPJ_UINT32)opj_int_min(dwt->sn, dwt->dn - a),

diff --git a/src/lib/openjp2/dwt.h b/src/lib/openjp2/dwt.h
index 89c859cb9727add166534fad4de568a52d7ae7ed..215061e6b9cf010da87b652b9a5f65f212e7f84b 100644 (file)
--- a/src/lib/openjp2/dwt.h
+++ b/src/lib/openjp2/dwt.h
@@ -73,12 +73,6 @@ OPJ_BOOL opj_dwt_decode(opj_tcd_t *p_tcd,
                         opj_tcd_tilecomp_t* tilec,
                         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.
 @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_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
 @param level Level of the wavelet function

diff --git a/src/lib/openjp2/t1.c b/src/lib/openjp2/t1.c
index 0787dce8df0980ebb693783c5139c9051fcebeca..937f420ad7e4846860bdf0c0b8e94613eafb5a07 100644 (file)
--- a/src/lib/openjp2/t1.c
+++ b/src/lib/openjp2/t1.c
@@ -1426,7 +1426,11 @@ static OPJ_FLOAT64 opj_t1_getwmsedec(
     if (qmfbid == 1) {
         w2 = opj_dwt_getnorm(level, orient);
     } else {    /* if (qmfbid == 0) */
+        const OPJ_INT32 log2_gain = (orient == 0) ? 0 :
+                                                (orient == 3) ? 2 : 1;
         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);

diff --git a/src/lib/openjp2/tcd.c b/src/lib/openjp2/tcd.c
index a84ec0639e2b38ba1f57af6c24a295120f2c202b..02fb11db23d764e08bfecd8aa69d283ca1b468ff 100644 (file)
--- a/src/lib/openjp2/tcd.c
+++ b/src/lib/openjp2/tcd.c
@@ -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_event_mgr_t* manager)
 {
-    OPJ_UINT32(*l_gain_ptr)(OPJ_UINT32) = 00;
     OPJ_UINT32 compno, resno, bandno, precno, cblkno;
     opj_tcp_t * l_tcp = 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 l_level_no;
     OPJ_UINT32 l_pdx, l_pdy;
-    OPJ_UINT32 l_gain;
     OPJ_INT32 l_x0b, l_y0b;
     OPJ_UINT32 l_tx0, l_ty0;
     /* 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_res = l_tilec->resolutions;
         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);*/
 
         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;
 
             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 );*/
 
                 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);
-                numbps = (OPJ_INT32)(l_image_comp->prec + l_gain);
+                {
+                    /* Table E-1 - Sub-band gains */
+                    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
                  * procedure" of the standard */
                 l_band->numbps = l_step_size->expn + (OPJ_INT32)l_tccp->numgbits -

diff --git a/tests/nonregression/test_suite.ctest.in b/tests/nonregression/test_suite.ctest.in
index 3cf4dfc2371463bdc1565b9ced04a2aaaf75ea40..32f87d37f036fadeffbc9a09ea84a0d9595d1916 100644 (file)
--- a/tests/nonregression/test_suite.ctest.in
+++ b/tests/nonregression/test_suite.ctest.in
@@ -32,16 +32,16 @@ opj_compress -i @INPUT_NR_PATH@/random-issue-0005.tif -o @TEMP_PATH@/random-issu
 # 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-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 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 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 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 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 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 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 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 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
 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:
 opj_compress -i @INPUT_NR_PATH@/Bretagne2.ppm -o @TEMP_PATH@/Bretagne2_5.j2k -c [64,64]
 # issue 316