]> Creatis software - gdcm.git/blob - src/gdcmopenjpeg/libopenjpeg/dwt.c
ENH: A user does not know (and should not) what we internally use for representing...
[gdcm.git] / src / gdcmopenjpeg / libopenjpeg / dwt.c
1 /*
2  * Copyright (c) 2001-2002, David Janssens
3  * Copyright (c) 2002-2004, Yannick Verschueren
4  * Copyright (c) 2002-2004, Communications and remote sensing Laboratory, Universite catholique de Louvain, Belgium
5  * Copyright (c) 2005, Reiner Wahler
6  * All rights reserved.
7  *
8  * Redistribution and use in source and binary forms, with or without
9  * modification, are permitted provided that the following conditions
10  * are met:
11  * 1. Redistributions of source code must retain the above copyright
12  *    notice, this list of conditions and the following disclaimer.
13  * 2. Redistributions in binary form must reproduce the above copyright
14  *    notice, this list of conditions and the following disclaimer in the
15  *    documentation and/or other materials provided with the distribution.
16  *
17  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS `AS IS'
18  * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
19  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
20  * ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
21  * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
22  * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
23  * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
24  * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
25  * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
26  * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
27  * POSSIBILITY OF SUCH DAMAGE.
28  */
29
30 /*
31  *  NOTE:
32  *  This is a modified version of the openjpeg dwt.c file.
33  *  Average speed improvement compared to the original file (measured on
34  *  my own machine, a P4 running at 3.0 GHz):
35  *  5x3 wavelets about 2 times faster
36  *  9x7 wavelets about 3 times faster
37  *  for both, encoding and decoding.
38  *
39  *  The better performance is caused by doing the 1-dimensional DWT
40  *  within a temporary buffer where the data can be accessed sequential
41  *  for both directions, horizontal and vertical. The 2d vertical DWT was
42  *  the major bottleneck in the former version.
43  *
44  *  I have also removed the "Add Patrick" part because it is not longer
45  *  needed.  
46  *
47  *  6/6/2005
48  *  -Ive (aka Reiner Wahler)
49  *  mail: ive@lilysoft.com
50  */
51
52
53 #include "dwt.h"
54 #include "int.h"
55 #include "fix.h"
56 #include "tcd.h"
57 #include <stdlib.h>
58
59 #define S(i) a[(i)*2]
60 #define D(i) a[(1+(i)*2)]
61 #define S_(i) ((i)<0?S(0):((i)>=sn?S(sn-1):S(i)))
62 #define D_(i) ((i)<0?D(0):((i)>=dn?D(dn-1):D(i)))
63 /* new */
64 #define SS_(i) ((i)<0?S(0):((i)>=dn?S(dn-1):S(i)))
65 #define DD_(i) ((i)<0?D(0):((i)>=sn?D(sn-1):D(i)))
66
67 /* <summary>                                                              */
68 /* This table contains the norms of the 5-3 wavelets for different bands. */
69 /* </summary>                                                             */
70 double dwt_norms[4][10] = {
71   {1.000, 1.500, 2.750, 5.375, 10.68, 21.34, 42.67, 85.33, 170.7, 341.3},
72   {1.038, 1.592, 2.919, 5.703, 11.33, 22.64, 45.25, 90.48, 180.9},
73   {1.038, 1.592, 2.919, 5.703, 11.33, 22.64, 45.25, 90.48, 180.9},
74   {.7186, .9218, 1.586, 3.043, 6.019, 12.01, 24.00, 47.97, 95.93}
75 };
76
77 /* <summary>                                                              */
78 /* This table contains the norms of the 9-7 wavelets for different bands. */
79 /* </summary>                                                             */
80 double dwt_norms_real[4][10] = {
81   {1.000, 1.965, 4.177, 8.403, 16.90, 33.84, 67.69, 135.3, 270.6, 540.9},
82   {2.022, 3.989, 8.355, 17.04, 34.27, 68.63, 137.3, 274.6, 549.0},
83   {2.022, 3.989, 8.355, 17.04, 34.27, 68.63, 137.3, 274.6, 549.0},
84   {2.080, 3.865, 8.307, 17.18, 34.71, 69.59, 139.3, 278.6, 557.2}
85 };
86
87
88 /* <summary>                          */
89 /* Forward lazy transform (horizontal).  */
90 /* </summary>                            */ 
91 void dwt_deinterleave_h(int *a, int *b, int dn, int sn, int cas) {
92     int i;
93     for (i=0; i<sn; i++) b[i]=a[2*i+cas];
94     for (i=0; i<dn; i++) b[sn+i]=a[(2*i+1-cas)];
95 }
96
97 /* <summary>                             */  
98 /* Forward lazy transform (vertical).    */
99 /* </summary>                            */ 
100 void dwt_deinterleave_v(int *a, int *b, int dn, int sn, int x, int cas) {
101     int i;
102     for (i=0; i<sn; i++) b[i*x]=a[2*i+cas];
103     for (i=0; i<dn; i++) b[(sn+i)*x]=a[(2*i+1-cas)];
104 }
105
106 /* <summary>                             */
107 /* Inverse lazy transform (horizontal).  */
108 /* </summary>                            */
109 void dwt_interleave_h(int *a, int *b, int dn, int sn, int cas) {
110     int i;
111 /*    for (i=0; i<sn; i++) b[2*i+cas]=a[i];*/
112 /*    for (i=0; i<dn; i++) b[2*i+1-cas]=a[(sn+i)];*/
113     int* ai;
114     int* bi;
115     ai=a;
116     bi=b+cas;
117     for (i=0; i<sn; i++) {
118       *bi = *ai;  bi+=2;  ai++;
119     }
120     ai=a+sn;
121     bi=b+1-cas;
122     for (i=0; i<dn; i++) {
123       *bi = *ai;  bi+=2;  ai++;
124     }
125 }
126
127 /* <summary>                             */  
128 /* Inverse lazy transform (vertical).    */
129 /* </summary>                            */ 
130 void dwt_interleave_v(int *a, int *b, int dn, int sn, int x, int cas) {
131     int i;
132 /*    for (i=0; i<sn; i++) b[2*i+cas]=a[i*x];*/
133 /*    for (i=0; i<dn; i++) b[2*i+1-cas]=a[(sn+i)*x];*/
134     int* ai;
135     int* bi;
136     ai=a;
137     bi=b+cas;
138     for (i=0; i<sn; i++) {
139       *bi = *ai;  bi+=2;  ai+=x;
140     }
141     ai=a+(sn*x);
142     bi=b+1-cas;
143     for (i=0; i<dn; i++) {
144       *bi = *ai;  bi+=2;  ai+=x;
145     }
146 }
147
148
149 /* <summary>                            */
150 /* Forward 5-3 wavelet tranform in 1-D. */
151 /* </summary>                           */
152 void dwt_encode_1(int *a, int dn, int sn, int cas)
153 {
154   int i;
155
156   if (!cas) {
157     if ((dn > 0) || (sn > 1)) {   /* NEW :  CASE ONE ELEMENT */
158       for (i = 0; i < dn; i++) D(i) -= (S_(i) + S_(i + 1)) >> 1;
159       for (i = 0; i < sn; i++) S(i) += (D_(i - 1) + D_(i) + 2) >> 2;
160     }
161   } else {
162     if (!sn && dn == 1)          /* NEW :  CASE ONE ELEMENT */
163       S(0) *= 2;
164     else {
165       for (i = 0; i < dn; i++) S(i) -= (DD_(i) + DD_(i - 1)) >> 1;
166       for (i = 0; i < sn; i++) D(i) += (SS_(i) + SS_(i + 1) + 2) >> 2;
167     }
168
169   }
170 }
171
172 /* <summary>                            */
173 /* Inverse 5-3 wavelet tranform in 1-D. */
174 /* </summary>                           */ 
175 void dwt_decode_1(int *a, int dn, int sn, int cas) 
176 {
177   int i;
178
179   if (!cas) {
180      if ((dn > 0) || (sn > 1)) { /* NEW :  CASE ONE ELEMENT */
181       for (i = 0; i < sn; i++) S(i) -= (D_(i - 1) + D_(i) + 2) >> 2;
182       for (i = 0; i < dn; i++) D(i) += (S_(i) + S_(i + 1)) >> 1;
183     }
184   } else {
185     if (!sn  && dn == 1)          /* NEW :  CASE ONE ELEMENT */
186       S(0) /= 2;
187     else {
188       for (i = 0; i < sn; i++) D(i) -= (SS_(i) + SS_(i + 1) + 2) >> 2;
189       for (i = 0; i < dn; i++) S(i) += (DD_(i) + DD_(i - 1)) >> 1;
190     }
191   }
192 }
193
194
195 /* <summary>                            */
196 /* Forward 5-3 wavelet tranform in 2-D. */
197 /* </summary>                           */
198 void dwt_encode(tcd_tilecomp_t * tilec)
199 {
200   int i, j, k;
201   int* a;
202   int* aj;
203   int* bj;
204   int w, l;
205
206   w = tilec->x1-tilec->x0;
207   l = tilec->numresolutions-1;
208   a = tilec->data;
209
210   for (i = 0; i < l; i++) {
211     int rw;         /* width of the resolution level computed                                                           */
212     int rh;         /* heigth of the resolution level computed                                                          */
213     int rw1;      /* width of the resolution level once lower than computed one                                       */
214     int rh1;      /* height of the resolution level once lower than computed one                                      */
215     int cas_col;   /* 0 = non inversion on horizontal filtering 1 = inversion between low-pass and high-pass filtering */
216     int cas_row;   /* 0 = non inversion on vertical filtering 1 = inversion between low-pass and high-pass filtering   */
217     int dn, sn;
218
219     rw = tilec->resolutions[l - i].x1 - tilec->resolutions[l - i].x0;
220     rh = tilec->resolutions[l - i].y1 - tilec->resolutions[l - i].y0;
221     rw1= tilec->resolutions[l - i - 1].x1 - tilec->resolutions[l - i - 1].x0;
222     rh1= tilec->resolutions[l - i - 1].y1 - tilec->resolutions[l - i - 1].y0;
223
224     cas_row = tilec->resolutions[l - i].x0 % 2;
225     cas_col = tilec->resolutions[l - i].y0 % 2;
226
227
228     sn = rh1;
229     dn = rh - rh1;
230     bj=(int*)malloc(rh*sizeof(int));
231     for (j=0; j<rw; j++) {
232       aj=a+j;
233       for (k=0; k<rh; k++)  bj[k]=aj[k*w];
234       dwt_encode_1(bj, dn, sn, cas_col);
235       dwt_deinterleave_v(bj, aj, dn, sn, w, cas_col);
236     }
237     free(bj);
238
239     sn = rw1;
240     dn = rw - rw1;
241     bj=(int*)malloc(rw*sizeof(int));
242     for (j=0; j<rh; j++) {
243       aj=a+j*w;
244       for (k=0; k<rw; k++)  bj[k]=aj[k];
245       dwt_encode_1(bj, dn, sn, cas_row);
246       dwt_deinterleave_h(bj, aj, dn, sn, cas_row);
247     }
248     free(bj);
249   }
250 }
251
252
253 /* <summary>                            */
254 /* Inverse 5-3 wavelet tranform in 2-D. */
255 /* </summary>                           */
256 void dwt_decode(tcd_tilecomp_t * tilec, int stop)
257 {
258   int i, j, k;
259   int* a;
260   int* aj;
261   int* bj;
262   int w, l;
263
264   w = tilec->x1-tilec->x0;
265   l = tilec->numresolutions-1;
266   a = tilec->data;
267
268   for (i = l - 1; i >= stop; i--) {
269     int rw;         /* width of the resolution level computed                                                           */
270     int rh;         /* heigth of the resolution level computed                                                          */
271     int rw1;      /* width of the resolution level once lower than computed one                                       */
272     int rh1;      /* height of the resolution level once lower than computed one                                      */
273     int cas_col;   /* 0 = non inversion on horizontal filtering 1 = inversion between low-pass and high-pass filtering */
274     int cas_row;   /* 0 = non inversion on vertical filtering 1 = inversion between low-pass and high-pass filtering   */
275     int dn, sn;
276
277     rw = tilec->resolutions[l - i].x1 - tilec->resolutions[l - i].x0;
278     rh = tilec->resolutions[l - i].y1 - tilec->resolutions[l - i].y0;
279     rw1= tilec->resolutions[l - i - 1].x1 - tilec->resolutions[l - i - 1].x0;
280     rh1= tilec->resolutions[l - i - 1].y1 - tilec->resolutions[l - i - 1].y0;
281
282     cas_row = tilec->resolutions[l - i].x0 % 2;
283     cas_col = tilec->resolutions[l - i].y0 % 2;
284
285     sn = rw1;
286     dn = rw - rw1;
287     bj=(int*)malloc(rw*sizeof(int));
288     for (j = 0; j < rh; j++) {
289       aj = a+j*w;
290       dwt_interleave_h(aj, bj, dn, sn, cas_row);
291       dwt_decode_1(bj, dn, sn, cas_row);
292       for (k = 0; k < rw; k++)  aj[k] = bj[k];
293     }
294     free(bj);
295
296     sn = rh1;
297     dn = rh - rh1;
298     bj=(int*)malloc(rh*sizeof(int));
299     for (j = 0; j < rw; j++) {
300       aj = a+j;
301       dwt_interleave_v(aj, bj, dn, sn, w, cas_col);
302       dwt_decode_1(bj, dn, sn, cas_col);
303       for (k = 0; k < rh; k++)  aj[k * w] = bj[k];
304     }
305     free(bj);
306
307   }
308 }
309
310
311 /* <summary>                          */
312 /* Get gain of 5-3 wavelet transform. */
313 /* </summary>                         */
314 int dwt_getgain(int orient)
315 {
316   if (orient == 0)
317     return 0;
318   if (orient == 1 || orient == 2)
319     return 1;
320   return 2;
321 }
322
323 /* <summary>                */
324 /* Get norm of 5-3 wavelet. */
325 /* </summary>               */
326 double dwt_getnorm(int level, int orient)
327 {
328   return dwt_norms[orient][level];
329 }
330
331 /* <summary>                             */
332 /* Forward 9-7 wavelet transform in 1-D. */
333 /* </summary>                            */
334 void dwt_encode_1_real(int *a, int dn, int sn, int cas)
335 {
336   int i;
337   if (!cas) {
338     if ((dn > 0) || (sn > 1)) {   /* NEW :  CASE ONE ELEMENT */
339       for (i = 0; i < dn; i++)
340        D(i) -= fix_mul(S_(i) + S_(i + 1), 12993);
341       for (i = 0; i < sn; i++)
342        S(i) -= fix_mul(D_(i - 1) + D_(i), 434);
343       for (i = 0; i < dn; i++)
344        D(i) += fix_mul(S_(i) + S_(i + 1), 7233);
345       for (i = 0; i < sn; i++)
346        S(i) += fix_mul(D_(i - 1) + D_(i), 3633);
347       for (i = 0; i < dn; i++)
348        D(i) = fix_mul(D(i), 5038);   /*5038 */
349       for (i = 0; i < sn; i++)
350        S(i) = fix_mul(S(i), 6659);   /*6660 */
351     }
352   } else {
353     if ((sn > 0) || (dn > 1)) {   /* NEW :  CASE ONE ELEMENT */
354       for (i = 0; i < dn; i++)
355        S(i) -= fix_mul(DD_(i) + DD_(i - 1), 12993);
356       for (i = 0; i < sn; i++)
357        D(i) -= fix_mul(SS_(i) + SS_(i + 1), 434);
358       for (i = 0; i < dn; i++)
359        S(i) += fix_mul(DD_(i) + DD_(i - 1), 7233);
360       for (i = 0; i < sn; i++)
361        D(i) += fix_mul(SS_(i) + SS_(i + 1), 3633);
362       for (i = 0; i < dn; i++)
363        S(i) = fix_mul(S(i), 5038);   /*5038 */
364       for (i = 0; i < sn; i++)
365        D(i) = fix_mul(D(i), 6659);   /*6660 */
366     }
367   }
368 }
369
370 /* <summary>                             */
371 /* Inverse 9-7 wavelet transform in 1-D. */
372 /* </summary>                            */
373 void dwt_decode_1_real(int *a, int dn, int sn, int cas)
374 {
375   int i;
376   if (!cas) {
377     if ((dn > 0) || (sn > 1)) {   /* NEW :  CASE ONE ELEMENT */
378       for (i = 0; i < sn; i++)
379        S(i) = fix_mul(S(i), 10078);   /* 10076 */
380       for (i = 0; i < dn; i++)
381        D(i) = fix_mul(D(i), 13318);   /* 13320 */
382       for (i = 0; i < sn; i++)
383        S(i) -= fix_mul(D_(i - 1) + D_(i), 3633);
384       for (i = 0; i < dn; i++)
385        D(i) -= fix_mul(S_(i) + S_(i + 1), 7233);
386       for (i = 0; i < sn; i++)
387        S(i) += fix_mul(D_(i - 1) + D_(i), 434);
388       for (i = 0; i < dn; i++)
389        D(i) += fix_mul(S_(i) + S_(i + 1), 12994); /* 12993 */
390     }
391   } else {
392     if ((sn > 0) || (dn > 1)) {   /* NEW :  CASE ONE ELEMENT */
393       for (i = 0; i < sn; i++)
394        D(i) = fix_mul(D(i), 10078);   /* 10076 */
395       for (i = 0; i < dn; i++)
396        S(i) = fix_mul(S(i), 13318);   /* 13320 */
397       for (i = 0; i < sn; i++)
398        D(i) -= fix_mul(SS_(i) + SS_(i + 1), 3633);
399       for (i = 0; i < dn; i++)
400        S(i) -= fix_mul(DD_(i) + DD_(i - 1), 7233);
401       for (i = 0; i < sn; i++)
402        D(i) += fix_mul(SS_(i) + SS_(i + 1), 434);
403       for (i = 0; i < dn; i++)
404        S(i) += fix_mul(DD_(i) + DD_(i - 1), 12994); /* 12993 */
405     }
406   }
407 }
408
409 /* <summary>                             */
410 /* Forward 9-7 wavelet transform in 2-D. */
411 /* </summary>                            */
412
413 void dwt_encode_real(tcd_tilecomp_t * tilec)
414 {
415   int i, j, k;
416   int* a;
417   int* aj;
418   int* bj;
419   int w, l;
420
421   w = tilec->x1-tilec->x0;
422   l = tilec->numresolutions-1;
423   a = tilec->data;
424
425   for (i = 0; i < l; i++) {
426     int rw;         /* width of the resolution level computed                                                     */
427     int rh;         /* heigth of the resolution level computed                                                    */
428     int rw1;      /* width of the resolution level once lower than computed one                                 */
429     int rh1;      /* height of the resolution level once lower than computed one                                */
430     int cas_col;   /* 0 = non inversion on horizontal filtering 1 = inversion between low-pass and high-pass filtering */
431     int cas_row;   /* 0 = non inversion on vertical filtering 1 = inversion between low-pass and high-pass filtering   */
432     int dn, sn;
433
434     rw = tilec->resolutions[l - i].x1 - tilec->resolutions[l - i].x0;
435     rh = tilec->resolutions[l - i].y1 - tilec->resolutions[l - i].y0;
436     rw1= tilec->resolutions[l - i - 1].x1 - tilec->resolutions[l - i - 1].x0;
437     rh1= tilec->resolutions[l - i - 1].y1 - tilec->resolutions[l - i - 1].y0;
438
439     cas_row = tilec->resolutions[l - i].x0 % 2;
440     cas_col = tilec->resolutions[l - i].y0 % 2;
441
442     sn = rh1;
443     dn = rh - rh1;
444     bj=(int*)malloc(rh*sizeof(int));
445     for (j = 0; j < rw; j++) {
446       aj = a + j;
447       for (k = 0; k < rh; k++)  bj[k] = aj[k*w];
448       dwt_encode_1_real(bj, dn, sn, cas_col);
449       dwt_deinterleave_v(bj, aj, dn, sn, w, cas_col);
450     }
451     free(bj);
452
453     sn = rw1;
454     dn = rw - rw1;
455     bj=(int*)malloc(rw*sizeof(int));
456     for (j = 0; j < rh; j++) {
457       aj = a + j * w;
458       for (k = 0; k < rw; k++)  bj[k] = aj[k];
459       dwt_encode_1_real(bj, dn, sn, cas_row);
460       dwt_deinterleave_h(bj, aj, dn, sn, cas_row);
461     }
462     free(bj);
463   }
464 }
465
466
467 /* <summary>                             */
468 /* Inverse 9-7 wavelet transform in 2-D. */
469 /* </summary>                            */
470 void dwt_decode_real(tcd_tilecomp_t * tilec, int stop)
471 {
472
473   int i, j, k;
474   int* a;
475   int* aj;
476   int* bj;
477   int w, l;
478
479   w = tilec->x1-tilec->x0;
480   l = tilec->numresolutions-1;
481   a = tilec->data;
482
483   for (i = l-1; i >= stop; i--) {
484     int rw;         /* width of the resolution level computed                       */
485     int rh;         /* heigth of the resolution level computed                      */
486     int rw1;      /* width of the resolution level once lower than computed one   */
487     int rh1;      /* height of the resolution level once lower than computed one  */
488     int cas_col;   /* 0 = non inversion on horizontal filtering 1 = inversion between low-pass and high-pass filtering */
489     int cas_row;   /* 0 = non inversion on vertical filtering 1 = inversion between low-pass and high-pass filtering   */
490     int dn, sn;
491
492     rw = tilec->resolutions[l - i].x1 - tilec->resolutions[l - i].x0;
493     rh = tilec->resolutions[l - i].y1 - tilec->resolutions[l - i].y0;
494     rw1= tilec->resolutions[l - i - 1].x1 - tilec->resolutions[l - i - 1].x0;
495     rh1= tilec->resolutions[l - i - 1].y1 - tilec->resolutions[l - i - 1].y0;
496
497     cas_col = tilec->resolutions[l - i].x0 % 2; /* 0 = non inversion on horizontal filtering 1 = inversion between low-pass and high-pass filtering */
498     cas_row = tilec->resolutions[l - i].y0 % 2; /* 0 = non inversion on vertical filtering 1 = inversion between low-pass and high-pass filtering   */
499
500     sn = rw1;
501     dn = rw-rw1;
502     bj = (int*)malloc(rw * sizeof(int));
503     for (j = 0; j < rh; j++) {
504       aj = a+j*w;
505       dwt_interleave_h(aj, bj, dn, sn, cas_col);
506       dwt_decode_1_real(bj, dn, sn, cas_col);
507       for (k = 0; k < rw; k++)  aj[k] = bj[k];
508     }
509     free(bj);
510
511     sn = rh1;
512     dn = rh-rh1;
513     bj = (int*)malloc(rh * sizeof(int));
514     for (j=0; j<rw; j++) {
515       aj = a+j;
516       dwt_interleave_v(aj, bj, dn, sn, w, cas_row);
517       dwt_decode_1_real(bj, dn, sn, cas_row);
518       for (k = 0; k < rh; k++)  aj[k * w] = bj[k];
519     }
520     free(bj);
521   }
522 }
523
524
525 /* <summary>                          */
526 /* Get gain of 9-7 wavelet transform. */
527 /* </summary>                         */
528 int dwt_getgain_real(int orient)
529 {
530   (void)orient;
531   return 0;
532 }
533
534 /* <summary>                */
535 /* Get norm of 9-7 wavelet. */
536 /* </summary>               */
537 double dwt_getnorm_real(int level, int orient)
538 {
539   return dwt_norms_real[orient][level];
540 }